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ABSTRACT 


This  thesis  analyzes  the  dynamic  stability  of  positively  buoyant  submersibles.  Six 
degree -of-frccdom  equations  of  motion  are  used  to  compute  steady  state  behavior  with 
motion  restricted  to  the  vertical  plane.  Steady  state  solutions  are  analyzed  for  various 
conditions  of  buoyancy  including  changes  in  (1)  the  amount  of  excess  buoyancy,  (2)  the 
location  of  the  center  of  buoyancy,  (3)  the  location  of  the  center  of  gravity,  as  well  as  (4) 
the  deflection  of  bow  and  stern  planes.  The  equations  of  motion  are  then  linearized  around 
these  steady  state  solutions  to  predict  dynamic  response  in  the  vertical  plane.  The  stability 
of  each  solution  is  then  determined  by  eigen  value  analysis.  The  study  then  expands  the 
analysis  to  include  all  six  degrees  of  freedom  (i.e.,  include  stability  analysis  in  the 
horizontal  plane).  Finally,  numerical  integration  methods  are  used  to  verify  the  results. 
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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

Controlling  emergency  ascent  situations  on  submersible  vehicles  such  as  dive  plane 
jam  recovery  is  of  concern  to  the  U.S.  Navy.  In  order  to  control  such  situations,  one  must 
first  be  able  to  predict  the  dynamic  response  of  positively  buoyant  submersibles. 

Dynamic  response  equations  of  motion  descibc  the  manuevering  characteristics  of 
submersible  vehicles  for  six  degrees  of  freedom.  These  equations  assume  constant 
coefficients  for  hydrodynamic  forces  and  moments  approximated  by  zero  frequeney  added 
mass  and  damping  terms  plus  the  quadratie  terms  for  drag  forees.  The  constant  coefficients 
vary  for  each  vehicle  and  are  dependent  on  such  things  as  vehicle  body  shape,  location  and 
magnitude  of  vehicle  weight,  location  and  magnitude  of  vehicle  buoyancy,  position  of  bow 
and  stern  planes,  position  of  rudder,  vehicle  speed,  vehicle  mass  characteristics,  vehicle 
hyd'-odynamic  coefficients,  propeller  rpm  and  control  surface  inputs.  This  thesis  uses  the 
equations  of  motion  and  hydrodynamic  coefficients  for  a  submerged  Mark  IX  Swimmer 
Delivery  Vehicle  (SDV)  developed  by  Smith,  Crane,  and  Summey  [Reference  1:11-16,21- 
3 1 J  to  forecast  the  dynamic  stability  of  a  submersible  in  a  positively  buoyant  condition. 

This  study  begins  by  using  the  six  equations  of  motion  to  compute  the  steady  state 
behavior  of  a  submersible  vehicle  with  motion  restricted  to  the  vertical  plane.  The  steady 
state  solutions  in  the  vertical  plane  are  calculated  for  various  conditions  of  buoyancy 
including  changes  in  the  amount  of  excess  buoyancy,  the  location  of  the  center  of 
buoyancy,  the  location  of  the  center  of  gravity,  as  well  as  the  deflection  of  bow  and  stern 
planes.  The  SDV's  equations  of  motion  are  then  linearized  around  these  steady  state 
solutions  to  predict  dynamic  response  motion  in  the  vertical  plane  for  the  various  conditions 


of  buoyancy.  Several  solutions  are  computed  and  the  stability  of  each  solution  is 
determined  by  eigen  va'ue  analysis.  The  thesis  then  expands  the  analysis  to  include  all  six 
degrees  of  freedom  (i.e.  include  stability  analysis  in  the  horizontal  plane).  Finally, 
numerical  integration  methods  arc  used  to  verify  the  results. 

B.  EQUATIONS  OF  MOTION 

The  six  degree  of  freedom  equations  of  motion  for  the  submersible  vehicle  shown  in 
Appendix  A  were  taken  from  Smith,  Crane,  and  Summey  [Reference  1:11-16]. 
Differentiation  with  respect  to  time  is  denoted  by  a  dot  over  the  variable;  e.g.  “  = 

These  equations  are  referenced  to  a  right-hand  orthogonal  axis  system  fixed  in  the  body 
(vehicle)  as  shown  in  Figure  1 .  Since  these  equations  are  in  reference  to  a  body  fixed  axis 
system,  the  Euler  angles  of  pitch  (0),roll  (<j)),  and  yaw  QP)  are  used  to  specify  orientation 
with  respect  to  the  inertial  reference  system.  The  rotation  sequence  for  q),  0  and  and  the 
Euler  angle  rates  for  <1>,  0  and  ^  shown  in  Appendix  B  were  taken  from  Smith,  Crane,  and 
Summey  [Reference  1:18-20].  Major  variables  and  parameters  as  defined  by  Smith, 
Crane,  and  Summey  [Reference  1:7- 10[  are  given  below  : 

I.  Dynamic  Variables 

u,v,w  -  Linear  velocity  components  of  vehicle  with  respect  to  origin  of 

body  axes  system  relative  to  fluid. 

p.q.r  -  Angular  velocity  components  of  vehicle  with  respect  to  body 

axes  system  relative  to  inertial  reference  system. 

X,Y.Z  -  Hydrodynamic  force  components  along  body  axes. 

K.M.N  -  Hydrodynamic  moment  components  along  body  axes. 
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2.  Mass  Distribution  Parameters 

m  -  Mass  of  the  flooded  vehicle,  including  the  mass  of  the 

entrained  lluid. 


W  -  Weight  of  the  flooded  vehicle,  including  the  weight  of  the 

entrained  lluid  (  =  g  m  ;  where  g  is  the  acceleration  of  gravity). 

V  -  Displacement  volume  of  the  vehicle. 


B  -  Buoyancy  force  acting  on  the  vehicle  (  =  p  g  V  ).  This  is 

independent  of  the  inertial  mass  distribution  of  the  submersible 
vehicle,  including  whether  or  not  it  is  Hooded. 


XQ,yQ,Zj^  -  Coordinates  of  the  CG  (center  of  gravity)  in  the  body  axis 

system  (Figure  I).  These  will  depend  on  the  mass 
distribution  of  the  vehicle,  including  the  mass  of  the  entrained 
fluid. 


Xg-y  b-^B  '  Coordinates  of  the  CB  (center  of  buoyancy)  in  the  body  axis 

system  (Figure  1).  These  are  independent  of  the  mass 
distribution  system,  but  may  vary  with  the  addition  or  removal 
of  external  appendages. 

I^.I  .1^  -  Moments  of  inertia  about  the  body  system  axes,  including  the 

entrained  fluid. 

^xy’^x7’^y/  ■  inertia  about  the  body  system  axes,  including  the 

entrained  fluid. 

3.  Remaining  Parameters 

p  -  Mass  density  of  fluid  medium 


1  -  Reference  length  used  to  nondimensionalizc  the  hydrodynamic 

coefficients. 
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b(x),  h(x) 


^nosc’^tail 


-  Width  and  height  of  vehicle  in  its  xy  and  xz  planes, 
respectively,  at  location  x  measured  in  the  body  axes  system 
(Figure  1).  These  quantities  are  required  in  the  integrals 
defining  the  crossflow  forces  and  moments  in  the  equations  of 
motion,  and  are  tabulated  within  the  Steady  State  Computer 
Program  (Appendix  C). 

-  Sternplane,  bowplane  and  rudder  deflection  angles  in  radians 
(Figure  1). 
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II.  SYSTEM  SOLUTIONS  IN  THE  VERTICAL  PLANE 


A.  GENERAL 

In  the  steady  state  condition,  the  submersible  will  have  reached  constant  linear  and 
angular  velocities.  Therefore,  the  body  fixed  linear  accelerations  (u  ,  v,  w)  and  the  body 
fixed  angular  accelerations  (p  ,  q  ,  f)  will  be  zero.  Similarly,  the  vehicle  will  have  reached 
a  constant  angle  of  pitch  (0)  making  its  derivatives  (0)  equal  to  zero.  Since  this  analysis  is 
restricted  to  steady  state  solutions  in  the  vertical  plane,  the  angle  of  roll  (<!))  and  its 
derivative  (<p)  will  be  zero  (in  Chapter  IV,  the  case  where  the  angle  of  roll  ((|))  is  180 
degrees  will  be  discussed  during  the  numerical  integration  analysis).  The  angle  of  yaw  (^) 
and  its  derivative  (^)  will  likewise  be  zero  due  to  the  vertical  plane  restriction.  It  should  be 
noted  that  had  this  analysis  not  been  restricted  to  the  vertical  plane,  steady  state  yaw  (^) 
would  not  necessarily  be  zero,  thereby  allowing  the  angular  velocities  (p,  q,  r)  to  be  non¬ 
zero.  However,  since  this  analysis  was  restricted  to  the  cases  where  <I),  0  ,  ^  ,  p,  q  and  r 
are  all  zero,  the  equations  of  motion  for  six  degrees  of  freedom  for  the  steady  state 
condition  can  be  reduced  to: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

(W  -  B)  sin  0  =  Xvv  v2  +  X^w  ^2  +  X^^^  uvbj  +  uw  (  X^^^  65  +  X^^j,  65), 

(  ^6s6s  ^6b6b  ^6r6r  ^r^)j  “^^prop 

•  Lateral  (Sway)  Equation  of  Motion: 

-  (W-B)  cos  0  sin<()=  Yy  uv  +  Yyw  vw  +  Y^j.  u2  bj- 

Coy  h(x)(v)2  +Cdz  b{x)(w)2  ^  ^  dx 


f^nose 

•'xtail 
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Normal  (Heave)  Equalion  of  Motion 
-  ( W-B)  cos  0  cos  <p  =  uw  +  Zvv  v2  +  u2 


J^^nose 
xiail 


Coy  h(x)  (v)2  +Cd^  b{x)  (w)2 


— dx 
Uci(x) 


Roll  Equation  of  Motion: 

-  (yG  W  -  yB  B)  cos  0  cos  <p  +  (zqW  -  zbB)  cos  0  sin  (p 
=  Ky  uv  +  vw^  +  u2  Kprop 

Pitch  Equation  of  Motion: 

(xG  W  -  XB  B)  cos  0  cos  <p  +  (zg  W  -  zb  B)  sin  0  = 

Mw  uw  +  Mvv  v2  +  u2  (m^j.  6s  +  M^^^  6b) 


f^nosv 

•'xtail 


CDy  h(x)  (v)2  +C£f2  b(x)  (w)2  .  X  dx 


Yaw  Equation  of  Motion: 

-  (xG  W  -  XB  B)  cos  0  sin  <t»  +  (yG  W  -  yB  B)  sin  0  = 

Ny  uv  +  Nyw  vw  +  u2  6r'  +  u2  Nprop 
r^nose 

-  I  Coy  h(x)  (v)2  +  Cdz  Nx) (w)2  ^  V  x  dx 

Jxtaii  ' 


B.  CONDITIONS 

1 .  Defining  Additional  Terms 
a.  Excess  Buoyancy,  bB 

Excess  buoyancy  is  defined  as  6B  =  B  -  W  where  B  is  the  submersible's 
total  buoyancy  and  W  is  the  submersible's  total  weight. 
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b.  Longitudunal  Center  of  Buoyancy,  xqb 

The  longitudinal  center  of  buoyancy  is  defined  as  =  xq  -  Xg  where 
xq  is  the  longitudinal  center  of  gravity  with  respect  to  the  body  fixed  axis  and  Xg  is  the 

longitudinal  center  of  buoyancy  with  respect  to  the  body  fixed  axis. 

c.  Vertical  Center  of  Buoyancy,  Zq^ 

The  vertical  center  of  buoyancy  is  defined  as  ZQg  =  Zq  -  Zg  where  zq  is 
the  longitudinal  center  of  gravity  with  respect  to  the  body  fixed  axis  and  Zg  is  the 
longitudinal  center  of  buoyancy  with  respect  to  the  body  fixed  axis  (zQg  is  assumed  to  be 

positive). 

2 .  Assumed  Conditions 

a.  Lateral  Centers  of  Gravity,  yQ  ,  and  Buoyancy,  yg 

The  lateral  center  of  gravity  and  the  center  of  buoyancy  are  assumed  to  be 
on  the  same  centerline  plane  (y^  =  =  0). 

b.  Propeller  Speed, n  (revolutions  per  minute) 

The  propeller  speed  is  assumed  to  be  zero  (n  =  0). 

c.  Propeller  Coefficients,  Kp^gp  and  Nprop 

From  Smith,  Crane,  and  Summey  [Reference  1:30],  the  propeller 
coefficients  are  zero  (Kp^p  =  Np^p  =  0). 

C.  REVISING  THE  EQUATIONS  OF  MOTION 

Using  the  term  for  venical  center  of  buoyancy  (zQg),  the  expression  (zqW  -  ZgB) 

may  be  written  as  (zQgW  -  Zg5B).  Similarly,  using  the  term  for  longitudinal  center  of 
buoyancy  (xQg)  the  expression  (xqW  -  XgB)  may  be  written  as  (xQgW  -  xg5B).  Also, 
the  term  u^Xpj-op  may  be  written  as: 

U^Xprop  =  “^CDo(Tl2-l}=u2CDo[(^^SMpdj2.i]  ^  CDoA2n2  -  Cdou2 
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where  A  is  a  constant. 


Since  the  shaft  speed  (n)  is  zero,  the  expression  may  be  further  reduced  to 
^^^prop  =  ■  Cdou^.  Substituting  these  expressions  plus  the  term  for  excess  buoyancy 

(6B)  and  the  assumed  conditions  revises  the  equations  of  motion  for  the  six  degree  of 
freedom  system  as  follows: 


Longitudinal  (Surge)  Equation  of  Motion; 

-6B  sin  e  =  Xvv  v2  +  X^vv  w2  +  X^^j.  uvbf  +  uw  (  X^^^.  6s  +  6h) 

(  ^6s6s  +  ^drdr  ^r^)_  *  ^DO 


Lateral  (Sway)  Equation  of  Motion: 

6B  cos  6  sin<|)=  Yy  uv  +  Yyw 
f^nosc 

-  CDyh(x)v2+CDi,b(x)w2’^^J^dx 

•'xtail 

Normal  (Heave)  Equation  of  Motion: 

6B  cos  0  cos  (p  =  Zn^,  uw  +  Zvv  ^2  +  u2  (Z^j.  6s+Z^j,  6b) 
f^nosc  _ 

-I  CDyh(x)v2+CD^b(x)w2  ;^j^dx 

/Xiuil  ' 


Roll  Equation  of  Motion; 

(zgbW  -  zb6B)  cos  0  sin  <f>=  Ky  uv  +  Ky  w  vw 


Pitch  Equation  of  Motion: 

(xGB  W-xb6B)  cos  0  cos  <p  +  (zGB  ^  -  ZB  ^b)  sin  0  = 
Mw  uw  +  Myy  v2  +  u2  (m^j.  6s  +  M^j,  6b) 

f^nosc 


Cdv  h(x)v2  +CDzb{x)w2 


1 1  j 


x  dx 


•  Yaw  Equation  of  Motion: 

(-XGB  W  +  XB  6B)  cos  0  sin  <p  =  uv  +  Nyw  vw  +  u2  6j-| 

Z’^nose . 

-  I  Coy  h(x)  v2  +  Cd2  b{x)  w2  ^  x  dx 

Jxtaii 

These  six  equations  only  have  five  unknowns:  u,v,w,6,  and  <)).  Therefore,  two  of 
these  equations  must  be  dependent  and  additional  conditions  are  required  in  order  to  make 
the  number  of  equations  equal  the  number  of  unknowns. 


D.  ADDITIONAL  CONDITIONS 

The  next  condition  to  be  applied  to  the  vehicle  is  that  the  rudder  will  remain 
centerlined,  that  is  6r  =  0.  Since  the  solutions  of  interest  are  those  in  which  the  vehicle 
remains  within  the  vertical  plane,  it  can  be  further  speeified  that  the  linear  velocity  in  the 
transverse  direction  (v)  equals  zero.  Recalling  that  the  angle  of  roll  (<]))  has  been  previously 
assumed  to  equal  zero,  the  trigonometric  functions  of  <|>  can  be  identified  as  sintp  equals  zero 
and  costp  equals  unity.  Substituting  these  quantities  back  into  the  equations  of  motion, 

three  of  the  six  equations  of  motion  (sway,  roll,  and  yaw)  yield  trivial  solutions.  In 
addition,  the  cross-flow  velocity  term  for  the  heave  and  pitch  equations  can  be 

reduced  to: 

Ucf(x)  =  (v  +xr)^  +  (w  -  xq)^  =  |w| 

since  v,  r,  and  q  are  zero.  Furthermore,  since  is  constant,  it  can  be  taken  outside  the 

integral.  Therefore,  the  three  remaining  equations  can  be  wriiten  as: 

•  Longitudinal  (Surge)  Equation  of  Motion: 

-6Bsin9.  +  X 

(  ^6s6s  +  ^6b5b  ’  ^do  u- 
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Normal  (Heave)  Equation  of  Motion: 


J^XfiosC 

b{x) 

Hail 


dx 


Pitch  Equation  of  Motion: 

(xGB  W-xb^b)  cos  0  +  (zGB  W  -  ZB  6B)  sin  6  = 

pnose 

u2  (m^^  6s  +  M^j,  6b)+  Cdz  b{x)  x 

■'xtail 


Mu.’  uw  + 


dx 


E.  VERTICAL  PLANE  EQUATIONS  OF  MOTION 


f^noK  *'novc 

By  defining  the  terms  A^  as  the  I  b(x)  dx  and  x^  as  ^  I  b(x)  x  dx,  the  three 

rX|;,,|  r.Kfail 

remaining  equations  defining  motion  in  the  vertical  plane  can  be  written  as: 


Longitudinal  (Surge)  Equation  of  Motion: 

-  6B  sin  0  =  Xww  w2+  uw  (  6s  +  6h) 


Normal  (Heave)  Equation  of  Motion: 

6B  cos  0  =  Z^r  uw  +  “2(z^6s+Zab<'b)-CD.  w  |w!  Av, 


•  Pitch  Equation  of  Motion: 

(xgb  '^•XB  ^b)  cos  0  +  (zgbW-zb6B)  sin  0  = 

Mw  uw  +  u2  (m^^  6s  +  M^j,  6b)+  Cdz  w  |w|  xaA^ 

F.  COMPUTER  PROGRAM  DEVELOPMENT 

Taking  these  three  equations  which  describe  motion  in  the  vertical  plane,  solving  the 
first  two  for  sine  0  and  cosine  0,  respectively;  and  dividing  all  three  through  by  u^  yields: 


n 


Longitudinal  (Surge)  Equation  of  Motion: 

,  Xww(«l2+(^)(X„6s%^Xw6bfb) 


u- 


6B 


(x 


6s6s^s  + 


•  Cdo 


Normal  (Heave)  Equation  of  Motion: 
cos 
u^ 


i  fs-^Zbb  Ob)  -  CD.  "}}'  KJ 

“  oB  U'^ 


•  Pitch  Equation  of  Motion: 

(xGB  W  -  XB  hB)  +  (zQB  W  -  zb  6B)  = 

Mw  "  +  (Mjs  h  +  6b)+  Cdz  ™  xaA.. 

Now,  defining  the  quantity  as  w',  and  substituting  w'  back  into  the  three 

equations: 


if  w  is  positive: 


•  Longitudinal  (Surge)  Equation  of  Motion: 

sin.e  ^  .  .1  (w')“+  (w')(  6s  +  6h) 

6B  +  (  6s^  +  6b^)  -  Cdo 

•  Normal  (Heave)  Equation  of  Motion: 

=  1  Zw(W)t(Zj,.  6s+Zj,9  6b)  -Cd.(W)2A, 

u2  5B 

•  Pitch  Equation  of  Motion: 

(’^GB  ^  ■  ’^B  +  (zGB  W  *  ZB  = 

u2  u2 

Mw  (W)  +  (M^j.  6s  +  M^^,  6b)  +  Cdz  (w')^  xaA^ 

however,  if  w  is  negative: 


Longitudinal  (Surge)  Equation  of  Motion; 
sine^.  1  Xww(W)^+ (w')(  X^^5  6s  +  X^^j^6b) 
u2  6B  +(x^^bs^s^+  X^bbb  *  ^DO 
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Normal  (Heave)  Equation  of  Motion: 

=  1  Zw(W)  +  (z^  6s-tZj|,  6h)  +CDz(wf  a.. 
u-  5B 


Pitch  Equation  of  Motion: 

(xgb  ^  +  (i^GB  ^  ■  ZB  ^b) 

u2  u2 

Mw  (w')  +  (%s  +  Mbb  ^b) '  Cdz  (w')2  xaAv 


•  »  V/  \£  .  .  'Ill 

Substituting  the  equations  lor  ^2  ^^2  ^be  pitch  equation  yields  the 

following  expressions: 

if  w  is  positive: 

(xGBW-xBbB)  *  Zw(w')  +  (Z^^  6s+Z^j,  6b)  -Cd/(w’)2 

6B 

■  (^OB  W  -  ZB  6B)  ‘  ’^“'7  ‘V  ^  = 

bB  +  (  6s  +  X^^5h  )  -  Cdo 

Mw  (w’)  +  (m^j,  6s  +  M^f,  6b)  +  Cdz  (w’f  xaAw 


and  if  w  is  negative: 

(xobW-xbOB)  Z,i.(W)  +  (z^  fts+Z^j,  6h)  +CDz(wf  A. 

Xww(w')2+(W)(x„j^6s  +  X„(,^6b) 
+ 


6B 

-  (zGB  W-zb^b)  ^ 
6B 


(^6s6s^s"+  ^6b6b^b^)’ 
Mw  (w')  +  (m^j,  6s  +  M^j^  6b)  -  Cd2  (w')2  xaA^ 


'“DO 


Rearranging  these  two  equations  to  get  them  into  the  form:  A(w')-  +  B(w')  +  C  =  0.  the 
expressions  become: 


if  w  is  positive: 

\aK  +  (zgB  ^  -  ZB  J-  Xww 

6B  {w’)2 

+  (xqb  W  -  XB  6B)  ^  Cd2  A^. 

6B 

Mw  +  (zGB  W  -  ZB  6B)  *  (  65  +  6b) 

+  6B  i(w') 

-  (xgb  W  -  XB  6b)  -i-Zw 

6B  j 

(zGB  W  -  ZB  6B)  M  6^2  -  C„„)  ' 

+  oB 

+  (Mbs  ^s  +  Mbb  ^b)-(xGB  W-xb6B)  —  (z^^.  6s+Z^j,  6b) 

6B 

and  if  w  is  negative: 

-  Cdz  xaAw  +  (zGB  W  -  ZB  6B)  ^  X^w  t 

6B  |(w’)2 

-  (xGB  ^  ■  XB  ^b)  ^  Cdz  Av  I 
6B  J 

Mw  +  (zGB  W  -  ZB  6B)  -i  (  X^bs  ^b)  ’ 

+  ;(w') 

-  (xGB  W  -  XB  ^b)  J-Zw  I 

6B  ' 

(zGB  W  -  ZB  6B|  i(  X^(,s  6^^  .  C„,) 

+  6B 

+  (Mbs  +  6b)-(xGBW-xB6B)  ^(z^^  6s+Z^b  ^b) 

6B 

These  quadratic  expressions  were  then  solved  using  the  equation: 

.  -B  ±  Tb^  -  4AC 


where: 


B  =6B  Mw  +  (zGB  W  -  ZB  6B)  (  X^bs  +  ^wbb  ^b) 

-(xgb  W-xb6B)  Zw 

c  =  (zGB  w  .  ZB  6B)  (  6,^  +  X;,^;,^  6b^  -  Cm,) 

+  6B(m^j,  6s  +  M^j,  bb)-(xGB  W-xb6B)  (z^^.  6s+Z^j^  6b) 


if  w  is  positive: 

A  =  6B(Cdz  xaA.v)  +  (zgb  W  -  zb  6B)  Xww  +  (xGB  -  xB  ^b)  Cdz  A^; 
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and  if  w  is  negative: 

A  =  -  6B{Cd^  xaA.)  +  (zgB  W  -  zg  6b)  -  (xcg  W  -  xg  6B) 

The  value  ol  B  was  determined  using  the  computed  values  of  w'  and  the  equation  for 


tangent  0: 


sin  0 

tan  0  = 

cos  0 


How'ever,  the  value  ot 


.COS0 

u- 


varied  depending  on  the  value  of  w  ',  w’hich  lead  to  two 


possible  solutions: 


tan  0  = 


Equation  for  tan  0  if  w'  is  Positive: 

-  ^w w  (w ' )-  -  (W  ’ )  (  6s  +  6h)  -  (  6s^  + 

Zw(w')++  (z^j,  6s+Z^h  ^b)  -CDz  Mw'Iw'I 


+  Cix» 


tan  0  = 


Equation  for  tan  0  if  w'  is  Negative: 

•  ■(“'■)(  +  ^w(ib '’!’]  ■  (  ^6s6s  +  ^bhilb  *!?5 

Z^,(w')  +  +  (Z^  6s+Z^I,  a^)  +Cd^  A>,<w'J|w'j 


+  Q)o 


In  either  case,  the  value  of  u^  was  computed  using  the  expression  derived  from  the 


surge  equation  of  motion: 


u^  = 


6B  sin  0 


-  X^^vv’  ( w ' )-  -  { W )  (  6s  +  6h)  -  (  X^j.5j.  6s^  +  6^^)  +  Coo 


This  leads  to  two  possible  solutions  for  u  (i.e.  u  =  ±  Vu^).  The  value  of  w  was 
computed  using  w  =  u  (w').  Combining  the  two  possible  solutions  of  u  with  the  two 
possible  values  for  w'  derived  from  the  quadratic  expression  lead  to  four  possible 
combinations  of  solutions  for  u  and  w.  The  computer  program  which  calculates  these  four 
possible  solutions  is  contained  in  Appendix  C.  It  is  an  interactive  program  designed  to 
allow  the  operaU' '  to  select  the  amount  of  excess  buoyancy  as  a  percentage  of  vehicle 
weight  (6B).  the  dellection  of  dive  planes  in  degrees  (6^,).  the  ratio  of  bow  planes  to  dive 
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weight  (6B),  the  dcneclion  of  dive  planes  in  degrees  (6^),  the  ratio  of  bow  planes  to  dive 
planes  the  location  of  and  Xg  from  body  fixed  axis  origin  as  a  percentage  of 

length,  and  the  location  of  z^^g  and  Zg  from  body  fixed  axis  origin  in  feet. 

G.  STEADY  STATE  RESULTS 

Figures  2,  3,  and  4  show  typical  steady  stale  solutions  for  surge  velocity,  heave 
velocity,  and  pitch  angle  as  a  function  of  dive  plane  angle.  The  two  cases  shown  vary  the 
location  of  the  longitudinal  center  of  buoyancy;  for  case  (a):  x^g  =  -  1  %  of  the  vehicle 

length  (L).  and  for  case  (b):  ''^GB  -  +  1  ^  L.  The  following  parameters  were  the  same  for 

both  cases;  excess  buoyancy  ,  6B  =  2  %  of  the  vehicle  weight  (W);  deflection  of  bow 
planes,  bj,  =  0;  location  of  horizontal  and  vertical  centers  of  buoyancy,  Xg  =  Zg  =  0;  and 
location  of  vertical  center  of  gravity,  z^g  =  0. 1  feel. 

All  runs  developed  four  solutions.  For  two  of  the  solutions,  the  magnitude  of  the 
surge  velocitities  were  large  while  the  magnitudes  of  the  associated  heave  velocities  were 
relatively  small.  This  has  been  descibed  by  Booth  |Ref  2:  297]  as  "predominantly  forward 
motion".  The  other  two  solutions  had  small  surge  velocity  magnitudes  and  large  heave 
velocity  magnitudes.  Booth  (Ref  3;  346]  referred  to  this  type  of  motion  as  "nearly  vertical 
ascents".  The  positive  or  negative  nature  of  the  velocities  arc  associated  with  the  value  of 
pitch  angle.  Positive  heave  velocities  are  associated  with  pitch  angles  greater  than  90 
degrees.  That  is  the  submersible  would  be  ascending  in  a  belly  up  orientation.  Although 
this  steady  state  analysis  computes  four  possible  solutions,  it  gives  no  indication  as  to 
which  of  the  solutions  are  stable  (if  any).  Dynamic  response  and  stability  criteria  will  be 
discussed  in  the  next  chapter. 
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Figure  2.  Steady  State  Vertical  Plane  Solutions  for  Surge 
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Figure  3.  Steady  State  Vertical  Plane  Solutions  for  Heave  Velocity  (w) 
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III.  DYNAMIC  STABILITY 


A.  GENERAL 

The  first  portion  of  this  chapter  uses  the  six  degree  of  freedom  equations  of  motion 
along  with  the  Euler  angle  rale  equations  for  the  derivatives  of  the  angles  of  pitch  and  roll 
(0  and  <[))  to  predict  the  dynamic  stability  when  movement  is  restricted  to  the  vertical  plane. 
These  equations  of  motion  are  then  linearized  around  the  vertical  plane  steady  state  nominal 
points  computed  in  chapter  II.  Eigen  value  analysis  is  then  used  to  compute  the  stability  of 
each  solution.  The  second  part  of  the  chapter  expands  this  analysis  to  include  all  six 
degrees  of  freedom  and  uses  the  same  steady  state  nominal  points  to  predict  the  dynamic 
stability  when  motion  is  not  restricted  to  the  vertical  plane. 

B.  RESTRICTING  MOTION  TO  THE  VERTICAL  PLANE 

Since  motion  is  restricted  to  the  vertical  plane,  the  body-fixed  transverse  veloctiy  (v) 
and  its  derivative  (v)  are  zero.  The  angles  of  roll  (<}))  and  yaw  ),  and  their  derivatives  (<t> 
and  are  also  zero.  We  will  continue  to  assume  that  the  lateral  center  of  gravity  and  the 
lateral  center  of  buoyancy  arc  on  the  same  centerline  plane  (yo  =  ye  =  0)>  and  the  rudder  is 
centerlincd  (6r  =0). 

C.  LINEARIZED  VERTICAL  PLANE  EQUATIONS  OF  MOTION 

Substituting  these  values  into  the  original  equations  of  motion  (Appendix  A),  yields 
trivial  results  for  three  of  the  six  equations:  Lateral  (sway),  Roll,  and  Yaw.  The  remaining 
equations  reduce  to  the  following  form: 
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Lxingitudinal  (Surge)  Equation  of  Motion: 

;  m  -  Xu  U  +  m  zq  q  =  i  -Coo  + 

+  +  X^5t,6b  uw  +  Xq^s^s  +  Xqftb^b!  uq+  [X^^^v'  w2 

+  Xvtq  -  m  wq+  "Xqq  +  mxQ  q2  -  (W  -  B)  sin  6 


Normal  (Heave)  Equation  of  Motion: 
m  -  w  -  m  xg  +  q  =  u^ 

+  uw  +  m  uq+  mzQ  q2 


f 


Cd/  b(x)(w-xq)-  ^  dx  +  (W  -  B)  cos  6 


Pitch  Equation  of  Motion: 

mz(;  u  -  Mw  +  m  XG  w  +  ly  -  Mq  q  =  M^^bs  +  M^^bb  u2+  M^^  uw 
+  Mq  -  mxG  uq-  m/c. 


JM  xnosc 

■  Ci>  b(x)(w-xq)2_  X  dx 

Ktail 


-  (xgW  -  xbB)  cos  0  -  (zcjW  -  zbB)  sin  0 


These  three  equations  which  describe  motion  restrieied  to  the  vertical  plane  are  functions  of 
four  variables  (u.w.q,0)  plus  their  dcriatives  (u,w,q,0).  The  equations  of  motion  and  the 
Euler  angle  rate  equation  for  0  were  linearized  using  the  following  generalized  procedure: 


BlI(iJ)u  +  Bll(i.2)w  +  Bll(i,3)q  +  B1  l(i,4)0  =  '  u  + 

:  OU  iv,  i  OW 


idfr  idl] 
w  +  .  q  + 


.:Wu 


I  dq 


0 


where  the  B1  l(i)'s  are  the  constants  associated  with  the  derivatives  of  the  variables,  the 
functions  f;  represents  the  right  hand  side  of  the  nonlinear  equations,  and  i  =  1  to  4 
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identifies  the  three  equations  of  motion  (surge,  heave,  and  yaw)  plus  the  Euler  angle  rate 
equation  for  0.  The  partial  derivatives  were  computed  as  follows: 


Partial  Derivatives  of  the  Longitudinal  (Surge)  EOM: 

-  ^  =  2(-Cdo  +  +  X^5b^b}wo  =  A1 1(1,1) 

all 

^  =  2X,^wo  +  (X,,;,,ds  +  X,.6b6b)uo  =  A1 1(1,2) 

=  (Xwq  -  m)wo+  (Xqfts^S  +  Xqftb^bjuo  =  A1 1(1,3) 

=(W-B)  cos  00  =  A1 1(1,4) 

a0 

Partial  Derivatives  of  the  Normal  (Heave)  EOM: 

=  Z^wq  +  2  (Zftj^ds  +  Z5t,6b)uo  =  A1 1(2,1) 

=  uoZv,  -  2Ci>AJw(^  =  A1 1(2,2) 

=  uo(Z^,  +  m)  +  2CnteA^^.XA!wo|  =  A1 1(2,3) 

=  -(W-B)sin  00  =  A1 1(2.4) 


af2 

au 

af2 

aw 

af2 

aq 

ai'z 

a0 


Note  on  differentiation  procedure:  The  cross-flow  velocity  term  (U^)  was 
reduced  to  Ut;f(x)  =  (v  +xr)^  +  (w  -  xq)^  ^  =  w^  =  |w  -  qj .  Allowing 
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the  integral  term  in  the  heave  equation  to  be  defined  as  Ij,  where: 
Ii  =  Q)/ 1  b(x)(w-xq)|w  -  xq!  dx 


f 


/  xtaii 

*\n 

,,  ,  =  2CJ 

d(w  -  q)  J 

1  b(x)0 

Mail 

^I, 

(  I 

1  ^(w  -  q)l 

dw 

\  d(w  -  q)  1 

1  dw  j 

dl,  J 

dll  \ 

1  ^(w  -q)| 

dq  1 

d(w  -  q)  1 

1  dq  f 

I 


J  XI 


2CdJw(,!  1  b(x)  dx  =  2Cdz  |wu1  A^. 


J*xnoNC 

X 

xui! 


b(x)  dx  =  -  2Cdz  |woI  xa  A„ 


Partial  Derivatives  of  the  Pitch  EOM: 

=  2(H,>  +  H,,6b)uo  +  MvvWo  =A1 1(3,1) 

r/u 

f  ^  u„  +2uuCdAvX>o1  =  A1 1(3,2) 

dw 

=  (H|  -  mx(i)uo  -(mz(,)w  o  -  2C|>/  Ia|wo|  =  A1 1(3,3) 

dq 

=  (xc,W  -  xijB)  sinHo  -  (zoW  -zbB)cos  Bq  =  A1 1(3.4) 

ae 


Nt)te  on  differentiation  procedure:  Again  the  cross-tlow  velocity  term  (U^.^)  was 
reduced  to  |w  -  q|.  Allowing  the  integral  term  in  the  pitch  equation  to  be  defined 
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as  where: 
I2  =  Cd/  I 


b(x)(w-xq)|w  -  xql  x  dx 


dl2 

c'Kw  -  q) 


b(x)  |w  -  xql  X  dx 


£  =(fl(w'-'q)  bWxdx  =  2C„,lw„!x^A. 

J  Kiatl 

Aq'  '(fl(w-'q)  K'^Sq'”)  =-2Cl>W)  b(x)  dx  =  -  2Cd,  |w„t  U 


Partial  Derivatives  of  the  Euler  Angle  Rate  Equation  for  6: 


=  0  =  An(4,l) 
du 


=0  =A1 1(4,2) 
dw  ^ 


^/‘’  =  1  =  A1 1(4,3) 


=0  =  All  (4,4) 

de 


The  constants  corresponding  to  the  derivatives  of  the  four  variables  (i.c.  the  left  hand 
side  of  the  four  equations)  arc  as  follows: 
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•  Constants  Associated  with  Derivatives  for  the  Longitudinal 

(Surge)  EOM: 

li  =i*  m  -  Xu  =  B1 1(1,1) 

q  =>  mzG  =  Bl  1(1,3) 

B1 1(1,2)  =  B1 1(1,4)  =  0 

•  Constants  Associated  with  Derivatives  for  the  Normal 

(Heave)  EOM: 

w  =>  m  -  Zw  =  B1 1(2,2) 

q  => -(mxo  +  =  B1 1(2,3) 

B1 1(2.1)  =  B1 1(2.4)  =  0 

•  Constants  Associated  with  Derivatives  for  the  Pitch  EOM: 

u  =>  mzo  =  B1 1(3,1) 

w  =5>  -  (Mvv  +  m  X(-,)  =  B1 1(3,2) 
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q 


ly  -  Mq  =  B1 1(3,3) 


B  11(3,4)  =  0 


Constants  Associated  with  Derivatives  for  the  Euler  Angle  Rate 
Equation  for  6  : 


e  =j>  1  =  B1 1(4,4) 

B1 1(4J)  =  B1 1(4,2)  =  B1 1(4,3)  =0 

These  expressions  can  be  arranged  in  a  matrix  format  to  form  the  linearized  equations 
of  motion  in  the  vertical  plane  about  the  nominal  steady  state  points.  The  matrix  format  is 
as  follows: 

Bll  X  X  1  =  All  X  XI 

where: 


All(l,l) 

A1 1(1,2) 

A1 1(1,3) 

All(l,4) 

All(2,l) 

A1 1(2,2) 

A1 1(2,3) 

A1 1(2,4) 

A1 1(3,1) 

A1 1(3,2) 

A1 1(3,3) 

A1 1(3,4) 

A1 1(4,1) 

A1 1(4,2) 

A1 1(4,3) 

A  11(4,4) 

B1I(1,1) 

Bll(l,2) 

Bll(l,3) 

Bll(l,4) 

Bll(2,l) 

B1 1(2,2) 

B 11(2,3) 

B1 1(2,4) 

B1 1(3,1) 

B 11(3,2) 

B 11(3,3) 

B1 1(3.4) 

Bll(4,l) 

B 11(4,2) 

B1 1(4,3) 

B1 1(4,4) 
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and 


u 

w 


XI  = 

q 

0 

D.  VERTICAL  PLANE  COMPUTER  PROGRAM  DEVELOPMENT 

The  matrix  format  of  the  linearized  equations  of  motion  was  used  in  the  computer 
program  shown  in  Appendix  D  to  predict  the  dynamic  stability  of  the  vehicle  with  motion 
restricted  to  the  vertical  plane.  The  program  is  interactive  in  that  it  allow's  the  operator  to 
select  which  of  the  four  data  files  from  the  steady  state  analysis  (Chapter  II)  will  be  used  to 
define  the  nominal  points  for  the  linearization  process.  An  eigen  system  subroutine  was 
used  to  find  eigen  values  and  eigen  vectors.  The  program's  output  was  called  the  degree  of 
stability  and  only  shows  the  largest  real  part  of  all  eigen  values.  The  stability  criteria  is 
such  that  the  degree  of  stability  must  be  negative  in  order  for  the  solution  to  be  stable. 

E.  VERTICAL  PLANE  DYNAMIC  RESULTS 

Of  the  four  possible  steady  state  solutions  computed  in  Chapter  II,  only  one  solution 
from  each  case  yielded  stable  characteristics.  There  were  some  cases  in  which  none  of  the 
solutions  were  stable  for  certain  ranges  of  parameters.  The  general  trend  of  the  linearized 
dynamic  results  are  fairly  well  demonstrated  by  the  two  cases  discussed  in  Chapter  II. 
Recalling  the  parameters  of  these  cases:  6B  =  2%,  ratio  of  bow  to  stern  planes  =  0, 

Zgb  =  0.1  feet,  and  Xg  =  Zg  =  0.  The  first  case  placed  the  longitudinal  center  of  gravity  aft 
of  the  longitudinal  center  of  buoyancy  (x^g  =  -  1  %),  and  the  second  case  placed  the 
longitudinal  center  of  gravity  forward  (x^g  =  +  I  %)•  Once  again  dive  plane  deflection 
angle  (6^)  was  varied  from  -  20  to  +  20  degrees  for  both  cases.  Figure  5  shows 
longitudinal  velocity  (u)  as  a  function  of  dive  plane  angle  (6^).  Case  One  (x^g  =  -  1  %) 
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showed  predominantly  forward  motion,  while  case  two  (x^g  =  -  1  %)  yielded  nearly 

vertical  ascents.  Figure  6  shows  vertical  motion  (w)  as  as  a  function  of  dive  plane  angle 
(6j).  The  results  concur  with  Figure  5,  case  one  shows  very  little  vertical  motion  while 

case  two  demonstrates  a  larger  value.  It  is  interesting  to  note  that  vertical  motion  (w)  for 
case  one  is  positive  for  dive  plane  angles  between  -  20  and  -4  degrees.  The  case  one  values 
of  0  shown  in  figure  7  for  dive  plane  angles  between  -  20  and  -4  degrees  concur  with  this 
observation.  The  values  of  0  greater  than  90  degrees  indicate  the  submersible  is  ascending 

in  the  belly  up  position.  The  stability  in  the  vertical  plane  is  shown  in  figure  8.  Degree  of 
stability  (e)  is  shown  as  a  function  of  dive  plane  angle  (6^).  This  figure  shows  both  cases 


i 

2f- 


-20  -15  -10  -5  0  5  10  15  20 

dalta  B 

Figure  5.  Stable  Vertical  Plane  Solutions  for  Surge  Velocity 
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theta 


-20  -15  -10  -5  0  5  10  15  20 


delta  s 

Figure  6.  Stable  Vertical  Plane  Solutions  for  Heave  Velocity 


-20  -15  -10  -5  0  5  10  15  20 

delta  ■ 

Figure  7.  Stable  Vertical  Plane  Solutions  for  Pitch  Angle  (0) 
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0 


Figure  8.  Degree  of  Stability  in  the  Vertical  Plane 


F.  LINE4RIZ^\TI0N  OF  FULL  EQUATIONS  OF  MOTION 

Referring  back  to  the  complete  equations  of  motions  shown  in  Appendix  A,  these  six 
equations  which  describe  motion  for  the  six  degree  of  freedom  system  are  functions  of 
eight  variables  (u,v,w,p.q,r,<p,6)  plus  their  deriatives  (u,v,w,p,q,f,<}),0).  The  six 
equations  of  motion  (Appendix  A)  along  with  the  Euler  angle  rates  (Appendix  B)  were 
linearized  as  folows: 

bjlu  +  bj2v  +  bj3w  +  bj4p  +  bj5q  +  bj6f  +  bjVcp  +  bj89  = 


^g. 


5J 

du  ju,, 


U  + 


V  + 


0 


dw 


w  + 


P  + 


dgi 


y-v  .  "P  Jpo  .  Jqo 


q  + 


L  dr  jfd 


r  + 


L  d4'  .<1,0 


.M-  e 
-  ''S  je. 


where  the  bj's  are  the  constants  associated  with  the  derivatives  of  the  variables,  the 
functions  gj  represents  the  right  hand  side  of  the  nonlinear  equations,  and  j  =  1  to  8 
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identifies  the  six  equations  of  motion  plus  the  Euler  angle  rate  equations  for  <p  and  0, 
respectively.  The  partial  derivatives  of  these  equations  were  computed  as  follows: 


Partial  Derivatives  of  Lxrngitudinai  (Surge)  EOM: 

' =(^vv6s  ^vvt)b  &h}wo  +2(x5sfts  6b  )uo  -  2Ci)ouo  =  al  1 

OU 

=  2X,^.W„  +  (x,^,,  6s  +  X,.5h  6b)u(,  =  al3 

(7W 

=  -mwo  +  XwqWo  +  (Xqf,s  6s  +  Xq5f,  6b)uo  =  al5 
=-(W-B)cose  =  al8 


=0  =  312  =  0  =  al4  -^5*-  =  0  =  alb 

dv  dp  dr 


*«‘-  =  0  =  al7 
d(j) 


Partial  Derivatives  of  Lateral  (Sway)  EOM: 

=YvU(i  +  YvwWo  -  CdAv  |woI  =  a22 
dv 

Tp'  ^ 

=  -  muo  +Yr  uo+  Ywr  Wo  -  Ci)zAv,x^lwo|=  a26 


=  (W  -  B)  COS0  =  a27 
dcp 
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^K=0  =  a21  *S^  =  0=a23  ^2=0=a25  -^-  =  0=a28 

f)u  dw  dq  cl0 


Note  on  differentiation  procedure:  The  cross-flow  velocity  term  (U^f)  is  given 

by:  Ucf(x)  =  (v  +xrp  +  (w  -  The  integral  term  in  the  sway  equation 

was  defined  as  where: 


COy  h(x)  (v+xr)2  +Cd^  b(x)  (w-xq)2  ,  dx 


f^nosc 

•'xtail 

J'^nose 

">(Scf)) 

xtail 


dx 


dh 

-1  ^'1 

*  1 

'y+jeru 

dv 

\  r»v  1 

'  Ucf  / 

dv  ' 

1  Ucf  / 

=  Ci>/Wo2  Aw  = 

=  Cdz  |wo|  Aw 

W(,^ 

=1 

V  +  xr  \  .  T 

'  /V 

'_+.xrUc 

dr 

\  (tr  H 

Ucf  /  ar\ 

Ucf  ) 

=  C|)z  A 

~  U|)/  AwXy\  1' 

W(,2 

dh 

=  ( 

1  V  +  xr) 

dvJ 

\  dw 

l\  Ucf  / 

^W  \  Ucf 

=1  '"M 

( ''  *  +i 

d  , 

V  +  xr)| 

^q 

Uq  I 

1  Ucf  ) 

aql 

Ucf  1 

Ucf-(v  +  xr)  ~-Ucf  t 
dv 

. ^ 


u 


cf  ’ 


'  dr  I 


=  0 
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Partial  Derivatives  of  Normal  (Heave)  EOM: 


dl4 

4  ^Ml 

w-  xq 

+  I 

w  -  xq\ 

dw 

1  (iw  II 

Ucf  1 

dw 

Ucf  } 

=  2  Ctte  A^w  -  xq)^-5  +  I  -  :  Ucf  -  (w  -  ; 

cf  L  J 

=  2  C02  AwWo^-L  +Cdz  AwWo2^^  jwol  -  Wo i  =  2  AJwqI 
Fol  Wo^L  FoIj 

because  -  =  ^Ucf^x)  =  (v  +xr)^  +  (w  -  xq)^j^'^2(w  -  xq)  = 


dl4 

dq 


U?) "  ‘  ( "isr) ' 

D.A.W„2’‘Al»i 


ai 


h\ 

Partial  Derivatives  of  Roll  EOM: 


dg4 

'6\ 


Kv  Uo  +  Kvw  Wo  =  a42 


=  -  m  zg  Wo  +  Kp  Uo  +  Kwp  wo  =  a44 

dp  ‘ 


dg4 

dr 


m  zg  Uo  +  Kr  Uo  +  Kwr  Wo  =  a46 


=  -  (  ZG  w  -  ZH  B)  COS0O  =  a47 
dip 


^g4 

du 


=  0=a41 


=  0  =  a43  =  0  =  a45 

dw  dq 


=  0  =  a48 

de 
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Partial  Derivatives  of  Pitch  EOM: 


-  =  Mw  Wo  +  2  6s  +  6b)uo  =  a51 
du 

=  My,  uo  +  2  Cdz  AwXaIwoI  =  a53 


dss  ■»  M 

-f  •  mXG  Uo  -  mZG  Wo  +  Mq  Uo 

dq 


-  2  Cdz  Ia|wo1  =  a55 


=  (xgW  -  xhB)  sin  0o  -  (zgW  -  zbB)  cos  0o  =  a58 
d0 


=  0  =  a52  =  0  =  a54  =  0  =  a56  =  0  =  a57 

dv  dp  dr  d(}) 


Note  on  differentiation  procedure:  The  integral  term  in  the  pitch  equation 
was  defined  as  where: 

Coy  h(x)  (v+xr)2  +Cdz  Mx)  (w-xq)2  j  ^ -  x  dx 


/■Xno; 

•'xtail 

'L  *"(  u„w) I 

^^tail  *'Xj3j| 


I4  X  dx 


dls  _  dl4 
dw  dw 

dis  _  dl4 

dq  dq 

dIs  _  dl4 
dv  dv 


xa=  2CdzAwXa|w(J 


Xa  =  Ci>z  AwXaVoI  =  Cdz  IaIwoI 


xa  =  0 


dls 

dr 


dl4 

dr 


xa=  0 
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Partial  Derivatives  of  Yaw  EOM: 


dv 


Nv  Uo  +  Nvw  Wo  -  Cdz  Aw  xaIwoI  =  a62 


dgo 

dp 


m  XG  Wo  +  Np  uo  +  Nwp  wo  =  a64 


dr 


m  XG  Uo  +  Nr  Uo  +  Nwr  Wo  -  Cdz  U  |wo|  =  a66 


=  (xg  w  -  XB  B)  COS0O  =  a67 
d<j) 


=  0  =  a6 1  =  0  =  a63  =  0  =  a65  =  0  =  a68 

c»u  dv,'  d(\  c»0 


Note  on  differentiation  procedure:  The  integral  term  in  the  yaw  equation 
;fined  ; 

^nose 


was  defined  as  where: 


l6  = 


Cdv  h(x)(v+xr)2  +CDzh{’t)(w-xq)2  ^  x  dx 


Jrxno! 

’^tail 

Jr^nosc  f^no; 

X., 


Ucttx) 


I3  x  dx 


dv 


dij 

dv 


Xa  =  Cdz  |wo|  Aw  Xa 


dr 


dr 


XA  =  Ci)^  AwXa^  |woI  =  Cdz  Ukol 


^w 


dh 

dw 


Xa  =  0 


li*  4'i  x^  =  o 

dq  dq 
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Partial  Derivatives  of  Euler  Angle  Rate  fur  <{>  : 


=  1  =  a74 

=  tan  6(1  =  a76 

dr 

=  0  =  a71 

=  0  =  a72 

-^2=0=  a73 

clu 

^v 

dw 

=  0  =  a75 

^  0  =  all 

-S7  =  0  =  a78 

dq 

d<]) 

Partial  Derivatives  of  Euler  Angle  Rate  for  6  : 


=  ,  =  a86 


=0  =  a81 

=  0  =  a82 

=  0  =  a83 

du 

d\ 

dw 

=  0  =  a86 

=  0  =  a87 

=  0  =  a88 

dr 

dtp 

ae 

=  0  =  a85 


The  constants  corresponding  to  the  derivatives  of  the  eight  variables  (i.e.  the  left  hand 
side  of  the  eight  equations)  are  as  follows: 
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*  Constants  Associated  with  Derivatives  for  the  LongitudinaJ 

(Surge)  EOM: 

li  =*•  (m-Xu)=bll  q  =>  (mzo)=bl5 

bl2  =  bl3  =  bl4  =  bl6  =  bl7  =  bl8=0 

*  Constants  Associated  with  Derivatives  for  the  Lateral  (Sway)  EOM: 

V  =>  (m  -  Y;)  =  b22  p  =>  (-mzo  -  Yp)  =  b24 

r  =>  (mxo  -  Yf)  =  b26 
b21  =  b23  =b25  =  b27  =  b28  =  0 

■  Constants  Associated  with  Derivatives  for  the  Normal  (Heave) 
EOM: 

w  =*.  (m  -  Zvi.)  =  b33  q  =>  (.  mxo  -  2^)  =  b35 

b3 1  =  b32  =  b34  =  b36  =  b37  =  b38 

*  Constants  Associated  with  Derivatives  for  the  Roll  COM: 

V  =>  (-mzG  -  Kv)  =  b42  r  =s>  (-Kr)=b46 

P  =>  (Ix-Kp)=b44 

b41  =  b43  =  b45  =  b47  =  b48 


Constants  Associated  with  Derivatives  for  the  Pitch  EOM: 

u  =>  (mz(,)  =  b5 1  w  =>  (-  mxo  -  M^)  =  b53 

q  =>  (Iy-Mq)=b55 
b52  =  b54  =  b56  =  b57  =  b58 


Constants  Associated  with  Derivatives  for  the  Yaw  EOM: 

V  =t>  (mxo  -  Nv)  =  b62  p  =>  ( -  Np)  =  b64 

f  =>  (I2-Nf)=b66 


b61  =  b63  =  b65  =  b67  =  b68  =  0 
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•  Constants  Associated  with  Derivatives  for  the  Euler  Angle  Rate 

Equation  for  (p  : 

<p  =>  1  =  b77 

b71  =  b72  =  b73  =  b74  =  b75  =  b76  =  b78  =  0 

*  Constants  Associated  with  the  Derivatives  for  the  Euler  Angle 

Rate  Equation  for  H  : 

H  =>  1  =  b88 

b81  =  b82  =b83  =  b83  =  b84  =  b85  =  b86  =  b87  =  0 

These  expressions  can  be  arranged  in  a  matrix  format  to  define  the  dynamic  equations 

of  motion  for  the  six  degree  of  freedom  system  linearized  about  the  nominal  steady  state 

points.  The  matrix  format  is  B  x  X  =A  x  X,  where: 

all  al2  al3  al4  al5  al6  al7  al8 

a2l  a22  a23  a24  a25  a26  a27  a28 

a31  a32  a33  a34  a35  a36  a37  a.38 

a41  a42  a43  a44  a45  a46  a47  a48 

A  = 

a51  a52  a53  a54  a55  a56  a57  a58 

a61  a62  a63  a64  a65  a66  a67  a68 

a71  a72  a73  a74  a75  a76  a77  a78 

a81  a82  a83  a84  a85  a86  a87  a88 


all 

0 

al3 

0 

al5 

0 

0 

al8 

0 

a22 

0 

a24 

0 

a26 

■all 

0 

a31 

0 

a33 

0 

a35 

0 

0 

a38 

0 

a42 

0 

a44 

0 

a46 

a47 

aO 

a51 

0 

a53 

0 

a55 

0 

0 

a58 

0 

a62 

0 

a64 

0 

a66 

a67 

0 

0 

0 

0 

1 

0 

a76 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 
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bll 

bl2 

bl3 

bl4 

bl5 

bl6 

bl7 

bl8 

b21 

b22 

b23 

b24 

b25 

b26 

b27 

b28 

b31 

b32 

b33 

b34 

b35 

b36 

b37 

b38 

b41 

b42 

b43 

b44 

b45 

b46 

b47 

b48 

b51 

b52 

b53 

b54 

b55 

b56 

b57 

b58 

b61 

b62 

b63 

b64 

b65 

b66 

b67 

b68 

b71 

b72 

b73 

b74 

b75 

b76 

b77 

b78 

b81 

b82 

b83 

b84 

b85 

b86 

b87 

b88 

-- 

bll 

0 

0 

0 

bl5 

0 

0 

0 

0 

b22 

0 

b24 

0 

b26 

0 

0 

0 

0 

b33 

0 

b35 

0 

0 

0 

0 

b42 

0 

b44 

0 

b46 

0 

0 

b51 

0 

b53 

0 

b55 

0 

0 

0 

0 

b62 

0 

b64 

0 

b66 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

xl  ~ 

u 

x2 

V 

x3 

w 

x4 

p 

x5  i  q 


x7  :  j  ^ 
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These  matrices  may  be  reordered  such  that  they  will  be  of  the  form: 


Bll  0  XI  All  0  X\  ] 
0  B22  1  ^  X2  J  ^  ^  ^22  j  1X2* 


This  is  accomplished  by  rewriting  the  X  matrix  such  that  XI  is  the  same  matix  used 
in  the  vertical  plane  analysis: 


X  = 


xl 

x3 

x5 

x8 

x4 

x7 

x2 

.x6„ 


u 

w 

q 

e 

p 

<i> 

V 


XI 
X2  ■ 


The  A  matrix  is  restructured  into  a  matrix  containing  four  4x4  matrices  with 
theA12  and  A21  matrices  containing  only  the  zero  element. 


all 

al3 

al5 

al8 

0 

0 

0 

0 

a31 

a33 

a35 

a38 

0 

0 

0 

0 

a51 

a53 

a55 

a58 
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0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

All 

0 

0 

0 

0 

0 

a44 

a47 

a42 

a46 

0 

A22 

0 

0 

0 

0 

1 

0 

0 

a76 

0 

0 

0 

0 

a24 

a27 

a22 

a26 

0 

0 

0 

0 

a64 

a67 

a62 

a66  _ 
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Similarly,  B  is  restructured  into  a  matrix  containing  four  4x4  matrices  with 


the  B12  and  B2 1  matrices  containing  only  the  zero  element. 


bll 

0 

bl5 

0 

0 

0 

0 

0  ' 

0 

b33 

b35 

0 

0 

0 

0 

0  1 

b51 

b53 

b55 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0  L 

Bll 

0  ; 

0 

0 

0 

0 

b44 

0 

b42 

b46 

1 

.  0 

B22 

0 

0 

0 

0 

0 

1 

0 

0  1 

0 

0 

0 

0 

b24 

0 

b22 

b26 

0 

0 

0 

0 

b64 

0 

b62 

b66  J 

As  discussed,  the  X 1  matrix  established  to  descibe  the  linearized  dynamic  stabilty  in  the 
vertical  plane  is  the  same  as  the  XI  matrix  within  X.  In  addition,  the  All  and  Bll 
matrices  from  the  vertical  plane  analysis  are  also  identical  to  those  in  X  .  That  is: 


All(l,l) 

A1 1(1,2) 

A1 1(1,3) 

A1 1(1,4)  ^ 

a44 

a47 

a42 

a46 

A1 1(2,1) 

A1 1(2,2) 

A1 1(2,3) 

A1 1(2,4)  i  ' 

1 

0 

0 

a76 

A1 1(3,1) 

A1 1(3,2) 

All  (3,3) 

A1 1(3,4) 

a24 

a27 

a22 

a26 

A1 1(4,1) 

A1 1(4,2) 

A  11(4,3) 

A1 1(4,4)  J  i 

a64 

a67 

a62 

a66  J 

Bll(l,l) 

Bll(l,2) 

B1 1(1,3) 

B1 1(1,4)  1  [ 

b44 

0 

b42 

b46  1 

Bll(2,l) 

B1 1(2,2) 

B  11(2,3) 

B1 1(2,4)  i 

0 

1 

0 

0  1 

B1 1(3,1) 

B1 1(3,2) 

B 11(3,3) 

B1 1(3,4)  "i 

b24 

0 

b22 

b26  ! 

Bll(4,l) 

B 11(4,2) 

B1 1(4,3) 

B1 1(4,4)  J  L 

b64 

0 

b62 

b66  .. 

TheA22.  B22,  and  X2  matrices  represent  the  additional  equations  necessary  to  describe 
the  linearized  motion  for  all  six  degrees  of  freedom;  henceforth  referred  to  as  the  horizontal 
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plane  contributions.  The  eigen  function  for  the  six  degree  of  freedom  model  is  computed 
by  taking  the  determinant  as  follows: 

All-XBll  0  !  ^0  =>(det  All-XBll':)(det  A22-XB22:)-0 

0  A22-AB22  : 

Since  the  eigen  function  will  be  the  product  of  these  two  determinants,  the 
resulting  eigen  values  will  merely  be  the  union  of  the  vertical  plane  eigen  values  and  the 
horizontal  plane  eigen  values.  The  significance  of  reducing  the  eigen  value  calculation  from 
an  8  by  8  matrix  problem  to  two  4  by  4  matrix  problems  is  not  in  the  computation  time 
saved.  But  rather  in  fact  that  now  the  horizontal  and  vertical  stabilities  have  been  separated 
and  identified. 

G.  COMPUTER  PROGRAM  DEVELOPMENT 

The  matrix  format  of  the  linearized  dynamic  response  equations  associated  with 
motion  in  the  horizontal  plane  was  added  to  the  computer  program  developed  previously 
(Appendix  D).  Once  again,  the  program  is  interactive  in  that  it  allows  the  operator  to  select 
which  of  the  four  data  files  from  the  steady  state  analysis  (Chapter  II)  will  be  used  to 
define  the  nominal  points  for  the  linearization  process.  An  eigen  system  subroutine  was 
used  to  find  eigen  values  and  eigen  vectors.  Two  outputs  were  added  to  the  program. 
First,  the  degree  of  stability  in  the  horizontal  plane,  and  next  the  degree  of  stability  of  both 
planes  (i.e.,  the  union  of  the  vertical  and  horizontal  degrees  of  stability).  Reminder: 
degree  of  stability  must  be  negative  in  order  for  the  solution  to  be  stable. 

H.  DYNAMIC  STABILITY  SOLUTIONS 

Qintinuing  with  the  same  cases  from  part  E,  the  horizontal  stability  for  the  two  cases 
are  shown  in  figure  9.  Case  two  (x^g  =  +  1  %)  is  stable  in  the  horizontal  plane  for  all 


values  of  dive  plane  angle  (dj;).  Whereas,  case  one  =  -  1  %)  is  unstable  in  the 
horizontal  plane  for  dive  plane  angle  (6^)  between  -20  and  -  9  degrees.  From  figure  7,  this 
corresponds  to  values  of  0  greater  than  140  degrees.  This  indicates  that  the  vehicle  is 
stable  (even  in  the  horizontal  plane)  for  values  of  0  greater  than  90  degrees.  A  submersible 
with  an  angle  of  pitch  greater  than  90  degrees  will  have  a  negative  metacentric  height  and 
will  therefore  be  statically  unstable.  However,  the  results  shown  in  figures  7  and  9  indicate 
that  the  vehicle  will  remain  dynamically  stable  is  such  a  condition.  This  'inverted 
pendulum'  type  stability  will  be  further  investigated  during  the  numerical  integration 
analysis  (Chapter  IV)  to  see  if  hydrodynamic  and  drag  forces  on  the  vehicle  can  actually 
cause  this  to  occur.  Figure  10  shows  the  combined  stabilities  for  the  horizontal  and  vertical 
planes.  It  should  be  noted  that  in  general  the  horizontal  plane  dictated  stabilty  for  the  cases 
considered. 
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Because  the  real  part  of  the  computed  eigen  values  must  be  negative  in  order  for  a 
solution  to  be  stable,  it  is  easy  to  identify  unstable  solutions.  However,  the  measure  of 
stability  for  a  stable  solution  is  harder  to  quantify.  The  magnitudes  shown  thus  far  for 
degree  of  stability  seem  very  small.  Does  this  mean  that  the  solutions  are  not  very  stable? 
In  an  attempt  to  answer  this  question,  the  values  for  degree  of  stability  for  a  neutrally 
buoyant  submersible  (6B  =  0  or  W  =  B)  were  computed.  Figure  11  shows  the  degree  of 
stability  in  (a)  the  horizontal  plane  and  (b)  the  vertical  plane  as  a  function  of  surge  velocity. 
This  figure  shows  that  the  magnitude  of  the  values  for  degree  of  stability  for  this  neutrally 
buoyant  case  are  indeed  of  the  same  order  of  magnitude  as  the  positively  buoyant  cases 
discussed  earlier. 


Sure:  Velocity  (u)  Surge  Velocity  (u) 

(a)  Horizontal  Plane  (b)  Vertical  Plane 


Figure  11.  Degree  of  Stability  as  a  Function  of  Surge  Velocity  (u)  for  a 

Neutrally  Buoyant  Submersible 
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IV.  SIMULATIONS  AND  DISCUSSION  OF  RESULTS 


A.  DYNAMIC  STABILITY  ANALYSIS  RESULTS 

The  dynamic  stability  analysis  included  in  this  Chapter  considers  stability  for  both 
planes  (i.e.  the  analysis  no  longer  breaks  stability  down  by  horizontal  or  vertical  plane). 
As  mentioned  in  Chapter  III,  horizontal  plane  stability  generally  dietates  the  stability  of  the 
vehicle. 

1.  Variations  in  Longitudinal  Center  of  Gravity  (x^g) 

Figures  12  through  15  show  how  ehanging  the  longitudinal  center  of  gravity 
(Xgg)  effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  amount  of 

excess  buoyancy  (6B)  is  two  percent  of  wieght  (W),  the  bow  plane  deflection  angle  (6|^)  is 
zero,  the  vertical  center  of  gravity  (i^Qg)  is  0.1  feet,  and  the  longitudinal  and  vertical  centers 
of  buoyancy  (Xgand  Zg)  are  zero.  The  longitudinal  center  of  gravity  (x^g)  is  varied  from 
-  2  to  +  2  percent  of  vehicle  length.  When  the  longitudinal  center  of  gravity  (x^g)  is  greater 

than  zero,  there  are  stable  solutions  for  the  full  range  of  dive  plane  angles.  However,  this 
is  not  true  when  the  longitudinal  center  of  gravity  (x^g)  is  less  than  zero.  The  range  of 

stable  solutions  is  restricted  when  XQg=  -  0.5  and  -  1 .0. 
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(b)  Xqb  =  0,  -  0.5,  -  1.0,  -  1.5  and  -  2.0  % 

Figure  12.  Stable  Surge  Velocity  (u)  Solutions  for  Variations  in  Xqb 
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(b)  xqb  =  0,  -  0.5,  .  1.0,  -  1.5  and  -  2.0  % 


Figure  15.  Degree  of  Stablity  for  Variations  in  x^g 
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2.  Variations  in  the  Amount  of  Excess  Buoyancy  (SB) 

Figures  16  through  19  show  how  changing  the  amount  of  excess  buoyancy 

(5B)  effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  bow  plane 
deflection  angle  (5^^)  is  zero,  the  longitudinal  center  of  gravity  (x^g)  is  -  0.5  percent  of 

vehicle  length,  the  vertical  center  of  gravity  (z^g)  is  0.1  feet,  and  the  longitudinal  and 
vertical  centers  of  buoyancy  (xg  and  Zg)  are  zero.  The  amount  of  excess  buoyancy  (SB)  is 

varied  from  0.5  to  2.9  percent  of  weight  (W).  The  only  case  where  there  are  stable 
solutions  for  the  full  range  of  dive  plane  angles  is  when  buoyancy  (SB)  is  0.5. 
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(b)  6B  =  2.5,  2.6,  2.7,  2.8,  2.9 

Figure  18.  Stable  Angle  of  Pitch  (6)  Solutions  for  Variations  in  6B 


3.  Variations  in  Vertical  Center  of  Gravity  (z^g) 

Figures  20  through  23  show  how  changing  the  vertical  center  of  gravity  (zQg) 

effects  the  d\  namic  response  of  the  submersible.  For  these  cases,  the  amount  of  excess 
buoyancy  (6B)  is  two  percent  of  weight  (W),  the  bow  plane  deflection  angle  (6^,)  is  zero, 
the  longitudinal  center  of  gravity  (x^g)  is  -  0.5  percent  of  vehicle  length,  and  the 
longitudinal  and  vertical  centers  of  buoyancy  (Xg  and  Zg)  are  zero.  The  vertical  center  of 
gravity  (z^g)  is  varied  from  0.05  to  0.30  feet.  The  general  trend  shown  in  these  figures  is 
that  the  larger  values  for  vertical  center  of  gravity  (z^g)  have  stable  solutions  over  a  larger 
range  of  dive  plane  angles. 

6 - - 

0.2S  ! 

5.5-  \  -i 


Figure  20.  Stable  Surge  Velocity  (u)  Solutions  for  Variations  in  ZQg 
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Figure  21.  Stable  Heate  Velocity  (w)  Solutions  for  Variations  in  2, 
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Figure  22.  Stable  Angle  of  Pitch  (0)  Solutions  for  Variations  in  z^. 
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Figure  23.  Degree  of  Stability  for  Variations  in 

4.  Variations  in  the  Longitudinal  Center  of  Buoyancy  (Xg) 

Figures  24  through  27  shou  how  changing  the  longitudinal  center  of  buoyancy 
(Xg)  effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  amount  of 
excess  buoyancy  (PB)  is  two  percent  of  weight  (W').  the  bow  plane  deflection  angle  (6^^)  is 
zero,  the  longitudinal  center  of  gras  ity  (x^g)  is  -  0.5  percent  of  vehicle  length,  the  vertical 
center  of  gras  ity  (zQg)  is  0.1  feet,  and  the  vertical  center  of  buoyancy  (zg)  is  zero.  The 
longitudinal  center  of  buoyancy  (Xgl  is  varied  from  -  9  to  +  2  percent  of  vehicle  length. 

The  general  trend  shown  in  these  figures  is  that  the  positive  values  for  longitudinal  center 
of  buoyancy  (Xg)  tend  to  have  stable  solutions  for  positive  dive  plane  angles.  On  the  other 

hand,  the  negative  values  for  longitudinal  center  of  buoyancy  (Xg)  tend  to  have  stable 
solutions  for  negative  dive  plane  angles. 
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5.  Variations  in  Bow  Planes  Deflection  Angle  (b^) 

Figures  28  through  3 1  show  how  a  non-zero  bow  planes  deflection  angle  (6^^) 

effects  the  dynamic  response  of  the  submersible.  For  these  cases,  the  amount  of  excess 
buoyancy  (6B)  is  two  percent  of  wieght  (W),  the  longitudinal  center  of  gravity  (x^g)  is 
-  0.5  percent  of  vehicle  length,  the  vertical  center  of  gravity  (z^g)  is  0.1  feet,  and  the 
longitudinal  and  vertical  centers  of  buoyancy  (xg  and  Zg)  are  zero.  The  bow  planes  are 

given  a  deflection  of  -  20  degrees.  The  significance  of  the  results  shown  in  these  figures  is 
that  for  certain  dive  plane  angles  (6^  =  -3  to  -12  degrees)  there  are  two  stable  solutions. 

11  - ! 


/ 

/ 


Figure  28.  Stable  Surge  Velocity  (u)  Solutions  for  a  Non-zero  Bow  Plane 

Angle  (6jj  =  -  20  degrees) 
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Figure  29.  Stable  Heave  Velocity  (w)  Solutions  for  a  Non»zero  Bow  Plane 

Angle  =  •  20  degrees) 
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Figure  30.  Stable  Pitch  Angle  (0)  Solutions  for  a  Non-zero  Bow  Plane 

•Angle  =  -  20  degrees) 
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Figure  31.  Degree  of  Stability  for  a  Non-zero  Bow  Plane  Angle 

=  -  20  degrees) 

B.  SIMULATIONS  USING  NUMERICAL  INTEGRATION  METHODS 

The  linearized  dynamic  response  results  were  verified  by  simulations  using  numerical 
integration  of  the  full  six  degrees  of  freedom  equations  of  motion  for  the  swimmer  delivery 
vehicle  (SDV).  Figure  32  shows  a  plot  of  angle  of  pitch  (0)  versus  time  for  the  center  of 
gravit}  forward  of  the  center  of  buoyancy  case  (xQg  =  +  1)  with  a  dive  plane  angle  (6^)  of 

-  15  degrees.  The  dotted  line  shows  the  linearized  results  from  figure  7,  the  solid  line 
shows  the  numerical  integration  results.  The  stead}  state  results  of  the  numerical 
integration  method  match  the  linearized  dynamic  results  exactly. 
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Figure  32.  Numerical  Integration  Solution  for  Angle  of  Pitch  (B)  when 
Center  of  Gravity  is  Forward  (x^.^  =  + 


Figure  .'>.1  shows  a  plot  ol'  angle  of  pitch  (B)  versus  time  I'or  the  center  of  gravity  alt 
o*'  the  center  ot  buoyancy  case  (x^g  =  *  I)  with  a  dive  plane  angle  (b^)  of  •  15  degrees. 

.■\ga'n  the  dotted  line  sho\ss  the  linearized  results  from  figure  7,  the  solid  line  shows  the 
numerical  integration  results.  And  once  again,  the  steady  state  results  of  the  numerical 
integration  method  match  the  linearized  dynamic  results  exactly.  However,  this  linearized 
dynamic  result  was  for  the  vertical  plane  only,  the  horizontal  plane  stability  analysis 
indicated  that  this  would  he  an  unstable  solution  (figure  9).  The  reason  for  this 
disagreement  in  the  results  is  investigated  by  adding  an  initial  angle  of  roll  to  the  numerical 
integration  analysis.  Adding  a  small  angle  of  initial  roll  =  1  degree)  caused  the  vehicle 

to  steady  out  at  137  degrees  vice  159  degrees  as  shown  in  figure  34.  This  initial  roll  angle 
also  caused  a  steadv  state  roll  anele  of  17  decrees  as  shown  in  figure  35.  In  turn,  this 
stcad\  state  roll  angle  caused  the  steadv  slate  yaw  velocity  (r)  shown  in  figure  36  (i.e. 


motion  is  nc'  longer  restricted  to  the  vertical  plane).  Figure  37  shows  a  plot  of  z  versus  x 
and  indicates  that  the  vehicle  is  taking  a  helical  ascent  as  dicussed  by  Booth  [Ref  2:  304- 
305].  Therefore,  the  numerical  integration  solution  resulting  in  a  steady  state  pitch  angle  of 
159  degrees  (figure  33)  is  unstable  in  the  horizontal  plane  as  predicted  by  the  linearized 
dynamic  response  analysis. 
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Figure  33.  Numerical  Integration  Solution  for  Angle  of  Pitch  (0)  when 
Center  of  Gravity  is  Aft  (x^g  =  -  19f) 
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Figure  34.  Numerical  Integration  Solution  for  Angle  of  Pitch  (6)  when 

and  Initial  Roll  Angle  (cp^)  is  1  degree 
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Figure  35.  Numerical  Integration  Solution  for  Angle  of  Roll  (<p)  when 
^GB  =  *  Initial  Roll  Angle  is  1  degree 
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Figure  36,  Numerical  Integration  Solution  for  Yaw  Velocity  (r)  when 

Initial  Roll  Angle  ((p^)  is  1  degree 
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Figure  37.  Numerical  Integration  Solution  for  Heave  Velocity  (z)  when 
Xqu  =  -  1^  and  Initial  Roll  Angle  is  1  degree 
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Some  question  still  remains  with  regards  to  the  measure  of  stabilty  of  the  'inverted 
pendulum'  solutions  predicted  by  both  the  linearized  dynamic  response  analysis  and  the 
numerical  integration  analysis.  The  linearized  dynamic  response  analysis  predicts  a  stable 
solution  for  the  case  when  center  of  gravity  is  placed  aft  of  center  of  buoyancy 
(Xgb  =  - 1  %)  and  the  dive  plane  angle  (6^)  is  -  7  degrees  (figure  10).  The  corresponding 

steady  state  value  for  pitch  angle  was  118  degrees  (figure  7).  A  random  persistent  roll 
disturbance  ((p^)  was  added  to  the  numerical  integration  model  and  the  results  are  shown  in 

figure  38.  The  solid  line  indicates  the  results  when  a  small  disturbance  is  added  ((j)^ 

centered  about  0. 1  degrees),  the  dashed  line  indicates  the  results  when  a  large  disturbance 
is  added  ((p^  centered  about  1 .0  degrees).  As  expected,  the  large  disturbance  caused  the 

vehicle  to  roll  over  as  shown  by  the  resulting  angles  of  roll  (<p)  in  figure  39.  However,  the 
vehicle  continued  to  remain  stable  during  small  constant  random  disturbances  in  the 
inverted  position  as  shown  by  the  resulting  angles  of  roll  (<p)  in  figure  40.  This  indicates 
that  indeed  these  'inverted  pendulum'  solutions  have  a  significant  measure  of  stability. 
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Figure  38.  Numerical  Integration  Solution  for  Angle  of  Pitch  (6)  when 
'gb  =  *  3  Persistent  Roll  Disturbance  is  Added 
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Figure  39.  Numerical  Integration  Solution  for  Angie  of  Roll  (cp)  when 
^CB  “  *  ®  Large  Persistent  Roll  Disturbance  is  Added 
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Figure  40.  Numerical  Integration  Solution  for  .\ngle  of  Roll  (cp)  when 
Xru  =  •  snd  a  Small  Persistent  Roll  Disturbance  is  Added 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS. 


The  steady  slate  analysis  resulted  in  four  possible  solutions  provided  that  vehicle 
motion  was  restricted  to  the  vertical  plane.  Analyzing  the  dynamic  stability  using  the 
steady  slate  results  as  nominal  points  generally  indicates  that  only  one  (if  any)  of  the 
four  possible  solutions  will  be  stable.  There  are  a  few  cases  where  two  solutions  are 
stable,  but  these  cases  (certain  non-zero  bow  plane  deflection  angles)  appear  to  be  the 
exception  and  not  the  rule. 

The  dynamic  stability  characteristics  of  submersibles  can  be  separated  with  respect  to 
vertical  plane  motions  (u,w,q,6)  and  horizontal  plane  motions  (v,p,r,(t)). 

It  is  possible  for  submersibles  to  be  dynamically  stable  with  respect  to  vertical  plane 
motion  in  the  inverted  (belly  up)  position  during  ascents  (’Inverted  Pendulum' 
stabilization). 

'Inverted  Pendulum'  stabilization  is  also  possible  in  the  horizontal  plane. 

Submersibles  are  able  to  maintain  this  inverted  orientation  (i.e.  ascend  belly  up 
without  rolling  over)  even  under  some  persistent  roll  excitation. 

As  a  recommendation,  the  dynamic  stability  analysis  should  be  expanded  to  include 
the  case  where  angle  of  roll  is  180  degrees,  and  the  case  where  angle  of  roll  is  non¬ 
zero  (i.e.  <p  neither  equals  zero  nor  180  degrees).  Analyzing  the  <p  =180  degrees  cases 
will  only  involve  changing  a  few  signs  with  regards  to  trigonometric  functions; 
however,  analyzing  the  non-zero  cases  will  require  significant  effort. 

Furthermore,  identifying  and  characterizing  different  stability  regions  over  a  range  of 
variations  of  the  system  parameters  should  be  the  matter  of  future  research. 
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APPENDIX  A 


SIX  DEGREE  OF  FREEDOM  EQUATIONS  OF  MOTION 

Source:  Smith,  Crane,  and  Summey  [Reference  1:11-16] 

1.  LONGITUDINAL  (SURGE)  EQUATION  OF  MOTION 

m  u  -  vr  -I-  wq  -  XQ  (q2  -I-  r^)  +  yQ  (pq  -  r)  +  zq  (pr  +  q)  = 

^pp  P“  +  ^qq  9?  +  ^rr  +  Xpj-  pr 
+  Xu  u  Xwq  wq  +  Xyp  vp  -t-  Xvr  vr 
+  uq  (  Xq5s  6s  +  Xq5b  6b)  +  Xfbr  urbf 

+  Xyv  ^2  -I-  Xv(/w  w*-  +  Xv6r  uv6i-  +  uw  (  ^wbs  6s  +  X^>5b  65) 

+  u*-  (  X5s5s  65“  +  Xbbbb  65^  +  X5j-5j-  6]-^)  -  (W  -  B)  sin  6 


2.  LATERAL  (S\NA^  )  EQUATION  OF  MOTION 

m  \  +  ur  -  wp  -t-  XG  (pq  +  r)  -  yo  (p“  +r2)  +  ZQ  “  P)  = 
Ypp  +  YjT  +  Ypq  pq+  Yqr  qr 
+  Yv  V  +  Ypup  -i-  Yr  ur  -f  Yvq  vq  -i-  Y^p  wp  +  Y^r  wr 
+  Yv  uv  +  Yvma  vw  -f-  Ybr  u2  bj- 

Coy  h(x)  (v-hxr)-  +Coz  b(x)  (w-xq)2  dx 

'  Ucftx) 

+  (W-B)  cos  6  sin({) 


^nose 


^lail 
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3.  NORMAL  (HEAVE)  EQUATION  OF  MOTION 

m  w  -  uq  +  vp  +  XQ  (pr  -  q)  +  Yq  (qr  +  p)  -  ZQ  (P^  +  q^I  = 
I  Zqq  +  Zpp  p-  +  Zpr  pr  +  Zjj 
+  Zw  w  +  Zquq  +  Zyp  vp  +  Zvr  vr  j 
+  uw  +  Zvv  V-  +  u2  (Z6s  6s+Z5b  6b)_' 

Coy  h(x)  (v+xr)2  +Cdz  b{x)(w-xq)2  j  dx 

Ucftx) 

+  (W-B)  cos  6  cos  (j) 


■^nose 


'xtail 


4.  ROLL  EQUATION  OF  MOTION 

^x  P  +  (Iz  “ly)  qr  +  ^xy  (pr  •  q)  •  lyz  (q^  ‘  r^) 

-  Ixz  (pq  +  r )  +  m  YG  (w  -  uq  +  vp)  -  zq  (v  +  ur  -  wp)  = 
Kp  p  +Kr  r  +  Kpq  pq  +  Kqr  qr 
+  Kv  ^  +  Kp  up  +  Kf  ur  +  Kyq  vq  +  Kwp  wp  +  Ky/r  wr 

+  Kv  uv  +  Kvw  +  (yG  ^  ■  YB  ®  ^ 

-  (zqW  -  zgB)  cos  0  sin  (j)  +  u2  Kpa,p 
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5. 


PITCH  EQUATION  OF  MOTION 

ly  ^  +  (Ix  ■  ly)  pr  -  ^xy  (<1^  +  p)  +  lyz  (PQ  -  r ) 


+  Ixz  (P”  -  r-)  -  m  XQ  (w  -  uq  +  vp)  -  zq  (u  -  vr  +  wp)  = 
Mq  q  +  Mpp  p-  +  Mpr  pr  +  Mfr 
+  MwW  +  Mq  uq  +  Myp  vp  +  Mvr  vr  : 

+  uw  +  Myv  v^  +  (M5s  65  +  M513  65)  . 


Xnose 


Coy  h(x)  (v+xr)2  +Cdz  b(x)  (w-xq) 


|2  i  (^-xq) 

■  Ucf(x) 


X  dx 


-  (xQ  W  -  xg  B)  cos  0  cos  (p  -  (zQ  W  -  zg  B)  sin  0 


6.  YAW  EQUATION  OF  MOTION 

^z  r  +  (ly  -  ^x)  Pq  -  Ixy  (p'  -  q')  -  lyz  (pr  +  q) 

+  Ixz  (qr  -  p)  -  m  XQ  (v  +  ur  -  wp)  -  yQ  (u  -  vr  +  wq)  = 

Np  p  +Nr  r  +  Npq  pq  +  Nqj  qr 
+  NyV  +  Np  up  +  Nf  ur  +  Nyq  vq  +Nvv'p  wp  +  Nwr  wr 
+  N'v  uv  +  Ny^y  vw  +  N^r  u-  6i- 

Coy  h(x)  (v+xr)2  +  Cdz  b(x)  (w-xq)2  x  dx 

Ucftx) 

-j-  (xQ  W  -  Xg  B)  cos  0  sin  <p  +  (yQ  W  -  yg  B)  sin  0  +  u-  Nprop 


75 


APPENDIX  B 


ROTATION  SEQUENCE  AND  EULER  ANGLE  RATES 

1 .  ROTATION  SEQUENCE  FOR  <p,  6  AND  ^ 

Smith.  Crane,  and  Summey  [Reference  1:8]  descibe  the  transition  from  body  fixed  to 
inertial  reference  frames  as  follows: 

Since  the  equations  of  motion  are  referred  to  an  axis  system  that  is  fixed  for  the 
SDV  (swimmer  delivery  vehicle),  and  thus  translates  and  rotates  with  it,  the 
orientation  and  position  of  the  moving  body  axis  system  relative  to  a  fixed  inertial 
reference  system  must  be  specified.  The  orientation  of  the  body  axis  system  with 
respect  to  the  inertial  reference  system  is  defined  by  the  standard  Euler  angles  ’1’ 
(yaw),  0  (pitch),  and  tp  (roll).  The  rotation  sequence  from  the  inertial  reference 

system  to  the  body  system  is  ,0  ,  and  tp  as  shown  in  Figure  B1  taken  from  Smith, 
Crane,  and  Summey  [Reference  1:18). 


2.  EULER  ANGLE  RATES  FOR  <p,0  AND  V 

The  Euler  angle  rates  used  along  with  the  six  equations  of  motion  (Appendix  A)  in 

order  to  completely  determine  the  motion  of  the  submersible  were  specified  by  Smith, 

Crane,  and  Summey  [Reference  1:20]  to  be: 

<p  =  p  +  q  sin  <p  tan  0  +  r  cos  <p  tan  0 

0  =  q  cos  <p  -  r  sin  <[) 
sin  <p  cos  o 

V  =  q -  + - 

cos  0  cos  0 
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(1)  Vehicle-Centered 
Gravity-Directed  System 
para  1 1 e 1  to  inertial 
axis  system. 


(2)  Horizontal  Flight 
Reference  System  derived 
from  XqYoZo  by  rotation  a- 
bout  Z  through  yaw  angle  rji. 


Reference  System  derived  (1*)  Vehicle  Body  Axis  Ref- 
from  X"Y"Z"  by  rotation  erence  System  derived  from 
apout  Y"  through  pitch  X’Y'Z'  by  rotation  about  X 


ang I e  9 . 


through  roll  angle 


Figure  Bl.  Unit  Sphere  Development  of  Euler  Angles 


onon  onon 


APPENDIX  C 


STEADY  STATE  COMPUTER  PROGRAM 


PROGRAM  STEADY 

STEADY  STATE  SOLUTIONS  IN  THE  VERTICAL  PLANE 

DIVE  PLANE  VARIATION 

REAL  L, MASS, IX, IY,IZ,IXZ,IYZ,IXY, LAMBDA 

REAL  KPDOT , KRDOT , KPQ , KQR , KVDOT , KP , KR , KVQ , KWP , KWR , KV 

REAL  KVW,KPN,KDB 

REAL  MQDOT , MPP , MPR , MRR , MWDOT , MQ , MVP , MVR , MW , MW , MDS 

REAL  MDB,NDRB 

REAL  NPDOT , NRDOT , NPQ , NQR , NVDOT , NP , NR , NVQ , NWP , NWR , NV 

REAL  NVW,NDRS 

DIMENSION  X(9) ,BR(9) ,HH(9) ,VEC1(9) 

GEOMETRIC  PROPERTIES  AND  HYDRODYNAMIC  COEFFICIENTS 

PI  =4.0*ATAN(1.0) 

WEIGHT=12000.0 
L  =  17.425 

RHO  =  1.94 

G  =  32.2 

CDO  =  0.0057 

MASS  =WEIGHT/G 
CDZ  =0.5*0.5*RHO 
XWW  =  1.710E-01*0.5*RHO*L**2 
XWDS  «  4 .600E-02*0.5*RHO*L**2 
XWDB  =  9.660E-03*0.5*RHO*L**2 
XDSDS  — 1.160E-02*0.5*RHO*L**2 
XDBDB  --8.070E-03*0.5*RHO*L**2 
CDO  =  CDO*0.5*RHO*L**2 

ZW  --3.020E-01*0.5*RHO*L**2 
ZDS  =-2.270E-02*0.5*RHO*L**2 
ZDB  — 2.270E-02*0.5*RHO*L**2 
MW  =  9.860E-02*0.5*RHO*L**3 
MDS  =-l .113E-02*0.5*RHO*L**3 
MDB  =  1 . 113E-02*0. 5*RH0*L**3 

OPEN( 21 , NAME-' STl .RES' , STATUS- ' NEW ' ) 
OPEN(22,NAME-'ST2.RES' , STATUS- ' NEW ' ) 

OPEN (23, NAME- ' ST3 . RES ' , STATUS- ' NEW ' ) 

OPEN( 24 , NAME-' ST4 .RES' , STATUS- ' NEW ' ) 

OPEN( 31 , NAME- 'COEF.DAT' , STATUS- ' HEW ' ) 
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o  n  o  n  no 


DEFINE  THE  LENGTH  X,  BREADTH  BR,  AND  HEIGHT  HH  TERMS 

X{1)  =  -105.9/12.0 

X(2)  =  -99.3/12.0 

X(3)  =  -87.3/12.0 

X(4)  =  -66.3/12.0 

X(5)  =  72.7/12.0 

X(6)  =  83.2/12.0 

X(7)  =  91.2/12.0 

X(8)  =  99.2/12.0 

X(9)  =  103.2/12.0 

BR(1)  «  0.00/12.0 

BR(2)  =  8.24/12.0 

BR(3)  =  19.76/12.0 

BR(4)  =  29.36/12.0 

BR(5)  =  31.85/12.0 

BR(6)  =  27.84/12.0 

BR(7)  =  21.44/12.0 

BR(8)  =  12.00/12.0 

BR(9)  =  0.00/12.0 

COMPUTE  AREA  AND  CENTROID 

CALL  TRAP( 9,BR,X,AREA) 

DO  9  1=1,9 
VEC1(I)=X(I)*BR(I) 

9  CONTINUE 

CALL  TRAP( 9,VEC1,X,XAA) 

XA=XAA/AREA 

WRITE  (*,1002) 

READ  (*,*)  DSMIND,DSMAXD, IDS 

DSMIN=DSMIND*PI/180 
DSMAX=DSMAXD*P 1/180 
WRITE  (*,1001) 

READ  (*,*)  RATIO 

WRITE  (*,1003) 

READ  ( *  ,  *  )  DELB 

DELB=DELB*WEIGHT/100 . 0 
WRITE  (*,1004) 

READ  ( *  ,  *  )  XGB 

XGB=XGB*L/100 . 0 
WRITE  (*,1005) 

READ  ( * , * )  ZGB 

WRITE  (*,1006) 

READ  ( * , * )  XB 

XB=XB*L/100 . 0 
WRITE  (*,1007) 

READ  (*,*)  ZB 

WRITE  (31,*)  RATIO,  DELB,  XGB,  ZGB,  XB ,  ZB 
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»> 


DO  1  1=1, IDS 

DS=DSMIN+ ( DSMAX-DSMIN ) * ( I-l )/( IDS-1 ) 

IF  (DELB. EQ. 0.0)  DELB«0 . 000001 
IF  (ZGB.EQ.0.0)  ZGB  -0.000001 
DB=RATIO*DS 

PX  -XGB*WEIGHT-XB*DELB 
PZ  -ZGB*WEIGHT-ZB*DELB 
DEN  -CDZ*AREA*(PX+XA*DELB) 

LAMBDA-MW*DELB-PX*ZW+PZ* ( XWDS+RATIO*XWDB ) *DS 
ALPHA  — PX* ( ZDS+RATIO*ZDB ) *DS-PZ*CD0+PZ* ( XDSDS+ 
RATIO*RATIO*XDBDB ) *DS**2+DELB* ( MDS+RATIO 
*MDB)*DS 
BETA  -PZ*XWW 
LAMB DA- LAMBDA/DEN 
ALPHA  -ALPHA  /DEN 
BETA  -BETA  /DEN 
C 
C 


A  - 

1 . 0+BETA 

B  - 

LAMBDA 

C  - 

ALPHA 

DET- 

B**2-4.0*A*C 

IF  (DET.LT.0.0)  GO  TO  2 
WP-(-B+SQRT(DET) )/(2.0*A) 
YY--XWW*WP**2-(XWDS*DS+XWDB*DB)*WP 
&  -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF  (WP.''E.0,0)  XX-ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
&  *WP*ABS(WP) 

IF  (WP.LT.0.0)  XX-ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
&  *WP*ABS(WP) 

THETA-ATAN2 ( YY , XX ) 

USQ-DELB *S IN ( THETA )/YY 
THETA-THETA* 18 0/P I 
DSD=DS*180/PI 
IF  (USQ.LT.0.0)  GO  TO  3 
IF  (WP.GE.0.0)  U-  SQRT(USQ) 

IF  (WP.LT.0.0)  U— SQRT(USQ) 

W-WP*U 

WRITE  (21,*)  DSD, THETA, U,W,WP 

3  WP=(-B-SQRT(DET) )/(2.0*A) 

YY=-XWW*WP**2- ( XWDS*DS+XWDB*DB ) *WP 
&  -(XDSDS*DS**2  +  XDBDB’"DB**2  )+CD0 

IF  (WP.GE.0.0)  XX-ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
&  *WP*ABS(WP) 

IF  (WP.LT.0.0)  XX-ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
&  *WP*ABS(WP) 

THETA-ATAN2 ( YY , XX ) 

USQ-DELB *S IN ( THETA )/YY 
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o  n 


DSD=DS*180/PI 
THETA-THETA*! 8 0/PI 
IF  (USQ.LT.O.O)  GO  TO  2 
IF  (WP.LT.0.0)  U»-SQRT(USQ) 

IF  (WP.GE.0.0)  U«  SQRT(USQ) 
W-WP*U 

WRITE  (22,*)  DSD, THETA, U,W,WP 


2  A  -  -1.0+BETA 

DET-  B**2-4.0*A*C 
IF  (DET.LT.0.0)  GO  TO  1 
WP«(-E+SQRT(DET) )/(2.0*A) 

yy— xww*wp*  *  2-  (  xwds  *ds-»-xwdb*db  )  *wp 

&  -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF  (WP.LT.0.0)  XX-ZW*WP+ZDS*DS+2DB*DB-CDZ*AREA 
*WP*ABS(WP) 

IF  (WP.GE.0.0)  XX-ZW*WP+ZDS*DS+ZDB*DB+CD2*AREA 
*WP*ABS(WP) 

THETA-ATAN2 ( yy , XX ) 

USQ-DELB*SIN(TKETA)/yy 
DSD-DS*180/PI 
THETA- THETA* 18 0/P I 
IF  (USQ.LT.O.O)  GO  TO  4 
IF  {WP.GE.0.0)  U— SQRT(USQ) 

IF  (WP.LT.0.0)  U-  SQRT(USQ) 

W-WP*U 

WRITE  (23,*)  DSD, THETA, U,W,WP 
4  WP«(-B-SQRT(DET) )/(2.0*A) 

yy--xww*wp**  2- ( xwds*ds+xwdb*db ) *wp 

&  -(XDSDS*DS**2+XDBDB*DB**2 )+CD0 

IF  (WP.LT.0.0)  XX-ZW*WP+ZDS*DS+ZDB*DB-CDZ*AREA 
*WP*ABS(WP) 

IF  (WP.GE.0.0)  XX-ZW*WP+ZDS*DS+ZDB*DB+CDZ*AREA 
*WP*ABS(WP) 

THETA-ATAK2 ( yy , XX ) 

usq-delb* SIN ( THETA )/yy 

DSD-DS*180/PI 
THETA-THETA* 18 0/PI 
IF  (USQ.LT.O.O)  GO  TO  1 
IF  (WP.GE.0.0)  U— SQRT(USQ) 

IF  (WP.LT.0.0)  U-  SQRT(USQ) 

W-WP*U 

WRITE  (24,*)  DSD, THETA, U,W,WP 


1 

CONTINUE 

STOP 

1001 

FORMAT  (  ' 

ENTER 

BOW 

PLANE  TO  DIVE  PLANE  RATIO' 

1002 

FORMAT  ( ' 

ENTER 

MIN, 

MAX,  AND  INCREMENTS  IN 

& 

DS  ( degrees } ' ) 

1003 

FORMAT  (  ' 

ENTER 

DELB 

(  %W  )  '  ) 

1004 

FORMAT  ( ' 

ENTER 

XGB 

(%L)') 

SI 


n  o  n  o 


1005  FORMAT  ('  ENTER  ZGB  (feet)') 

1006  FORMAT  ('  ENTER  XB  (%L)') 

1007  FORMAT  ('  ENTER  ZB  (feet)') 

END 

C 

SUBROUTINE  TRAP ( N , A , B , OUT ) 

NUMERICAL  INTEGRATION  ROUTINE  USING 
THE  TRAPEZOIDAL  RULE 

DIMENSION  A(1 ) ,B( 1 ) 

Nl-N-1 
OUT-0 . 0 
DO  1  I-1,N1 

OUTl-0 .5*(A(I)+A(I+1))*(B(I+1)-B(I)) 
OUT  -OUT+OUTl 
1  CONTINUE 
RETURN 
END 
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APPENDIX  D 


LINEARIZED  DYNAMIC  STABILITY  COMPUTER  PROGRAM 


C  PROGRAM  LINEARIZED  DYNAMIC  STABILITY 

C  10  20  30  40  50 

C2345678901234567890123456789012345678901234567890123456 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

DOUBLE  PRECISION  L , MASS , lA , IX , I Y , I Z 

DOUBLE  PRECISION  KPDOT , KRDOT , KPQ , KQR , KVDOT , KP , KR , 

&  KVO,KWP,KWR,KV, 

&  KVW , KPN , KDB , MQDOT , MPP , MPR , MRR , MWDOT , MQ , MVP , MVR , 

&  MW,MW,MDS,MDB, 

&  NDRB , NPDOT , NRDOT , NPQ , NQR , NVDOT , NP , NR , NVQ , NWP , NWR , 

&  NV,NVW,NDRS 
C 

DIMENSION  A1 ( 4 , 4 ) ,B1 (4 ,4 ) ,BETAl( 4 ) ,ALFRl ( 4 ) ,ALFI1 ( 4 ) 
DIMENSION  BB1(4,4) ,BB2(4,4),ZZZH4,4) ,ZZZ2(4,4) 
DIMENSION  A2{4,4) ,B2(4,4) ,BETA2(4) ,ALFR2(4) ,ALFI2(4) 
DIMENSION  WR1(4) ,WR2 ( 4 ) ,Wll ( 4 ) ,WI2(4) 

DIMENSION  X(9) ,BR(9) ,VEC1(9) ,VEC2(9) 

C 

C  GEOMETRIC  PROPERTIES 

C 

PI  =  4 .D0*DATAN(l.DO) 

WEIGHT*  12000.0 
IX  *  1760.0 

lY  *  9450.0 

IZ  =10700.0 

L  =  17.425 

RHO  =  1.94 

G  =  32.2 

CDO  =  0.0057 

MASS  =  WEIGHT/G 

CDZ  =  0.5*0.5*RHO 

CDO  =  CD0*0 . 5*RH0*L**2 

C 

C  SURGE  HYDRODYNAMIC  COEFFICIENTS 

C 

XPP  =  7 . 030E-03*0 . 5*RH0>^L**4 

XQQ  =-l . 470E-02*0 . 5*RH0*L**4 

XRR  =  4 . 010E-03*0 . 5*RH0*L*’*-4 

XPR  =  7 . 640E-04*0 . 5*RH0*L**4 

XUDOT  =-7 . 580E-03*0 . 5*RH0*L**3 

XWQ  =-l . 920E-01*0 . 5*RH0*L**3 


onn  nno  non 


XVP  =-3 . 240E-03*0 . 5*RH0*L**3 
XVE  =  1 . e90E-02*0 . 5*RHO*L**3 
XQDS  =  2 . 610E-02*0 . 5*RHO*L**3 
XQDB  =-2 . 600E-03*0 . 5*RHO*L**3 
XRDR  =-8 . 180E-04*0 . 5*RH0*L**3 
XW  =  5.290E-02*0.5*RHO*L**2 
XWW  -  1.710E-01*0.5*RHO*L**2 
XVDR  «=  1 .730E-03*0.5*RHO*L**2 
XWDS  «  4.600E-02*0.5*RHO*L**2 
XWDB  =  9 . 660E-03*0 . 5*RHO*L**2 
XDSDS  — 1 .160E-02*0.5*RHO*L**2 
XDBDB  *-8 . 070E-03*0 . 5*RHO*L**2 
XDRDR  =-l . 010E-02*0 . 5*RHO*L**2 
XRES  »  CD0*0 . 5*RHO*L**2 

SWAY  HYDRODYNAMIC  COEFFICIENTS 

YPDOT  =  1 . 270E-04*0 . 5*RHO*L**4 
YRDOT  =  1.240E-03*0.5*RHO*L**4 
YPC  =  4.125E-03*0.5*RHO*L**4 
YQR  --6 . 510E-03*0 . 5*RHO*L**4 
YVDOT  =-5 . 550E-02*0 . 5*RHO*L**3 
YP  =  3 . 055E-03*0 . 5*RHO*L**3 
YR  =  2 . 970E-02*0 . 5*RHO*L**3 
YVQ  -  2 . 360E-02*0 . 5*RHO*L**3 
YWP  =  2 . 350E-01*0 . 5*RHO*L**3 
YWR  --1 .880E-02*0.5*RHO*L**3 
YV  =-9 . 310E-02*0 . 5*RHO*L**2 
YVW  =  6 . 840E-02*0 . 5*RHO*L**2 
YDRS  -+2.270E-02*0.5*RHO*L**2 
YDRB  -+2 . 270E-02*0 . 5*RHO*L**2 

HEAVE  HYDRODYNAMIC  COEFFICIENTS 

ZQDOT  =-6 . 810E-03*0 . 5*RHO*L**4 
ZPP  =  1 .270E-04*0 . 5*RHO*L**4 
ZPR  =  6 . 670E-03*0 . 5*RHO*L**4 
ZRR  =-7.350E-03*0.5*RHO*L**4 
ZWDOT  =-2 . 430E-01*0 . 5*RHO*L**3 
ZQ  •=-!  .  350E-01*0 . 5*RHO*L**3 
ZVP  «-4 . 810E-02*0 . 5*RHO*L**3 
ZVR  -  4 . 550E-02*C . 5*RHO*L**3 
ZW  =-3 . 020E-01*0 . 5*RHO*L**2 
ZW  =-6 . 840E-02*0 . 5*RHO*L**2 
ZDS  =-2 . 270E-02*0 . 5*RH0*L*’"2 
ZDB  =-2 . 270E-02*0 . 5*RH0*L**: 

ROLL  HYDRODYNAMIC  COEFFICIENTS 

KPDOT  =-l . 010E-03*0 . 5*RH0*L**5 
KRDOT  =-3 . 370E-05*0 . 5*RHO*L**5 
KPQ  *-6 . 930E-05*0 . 5*RH0’^L**5 
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KQR  =  1.680E- 
KVDOT  =  1.270E- 
KP  =-1.100E- 
KR  =-8.410E- 
KVQ  =-5.115E- 
KWP  =-1.270E- 
KWR  =  1.390E- 
KV  -  3.055E- 
KVW  — 1.870E- 


02*0.5*RHO*L**5 
04*0 . 5*RHO*L**4 
02*0 . 5*RHO*L**4 
04*0 . 5*RHO*L**4 
03*0 . 5*RHO*L**4 
04*0 . 5*RHO*L**4 
02*0.5*RHO*L**4 
03*0.5*RHO*L**3 
01*0.5*RHO*L**3 


PITCH  HYDRODYNAMIC  COEFFICIENT 
MQDOT  »-l .680E-02*0.5*RHO*L**5 


MPP 

B 

5.260E- 

MPR 

B 

5.040E- 

MRR 

B. 

-2.860E- 

MWDOT 

B. 

-6.810E- 

MQ 

B. 

-6.860E- 

MVP 

B 

1 .180E- 

MVR 

S 

1.730E- 

MW 

= 

9.860E- 

MW 

B. 

-2 . 510E- 

MDS 

B 

-1 .113E- 

MDB 

S 

1.113E- 

05*0.5*RHO*L**5 
03*0 . 5*RHO*L**5 
03*0.5*RHO*L**5 
02*0 . 5*RHO*L**4 
02*0.5*RHO*L**4 
03*0.5*RHO*L**4 
02*0.5*RHO*L**4 
02*0.5*RHO*L**3 
02*0 . 5*RHO*L**3 
02*0 . 5*RHO*L**3 
02*0.5*RHO*L**3 


YAW  HYDRODYNAMIC  COEFFICIENTS 


NPDOT 

B-  3 

NRDOT 

=-3 

NPQ 

=— 2 

NQR 

=  2 

NVDOT 

=  1 

NP 

=  -8 

NR 

=  -l 

NVQ 

=  -9 

NWP 

=-l 

NWR 

=  7 

NV 

=  -7 

NVW 

=  -2 

NDRS 

=  -l 

NDRB 

=+l 

370E-05*0. 
400E-03*0. 
110E-02*0. 
750E-03*0. 
240E-03*0. 
405E-04*0. 
640E-02*0  . 
990E-03*0. 
750E-02*0. 
350E-03*0. 
420E-03*0. 
670E-02*0. 
113E-02*0. 
113E-02*0. 


5*RHO*L**5 

5*RHO*L**5 

5*RHO*L**5 

5*RHO*L**5 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**3 


DEFINE  THE  LENGTH  X  AND  BREADTH  BR  TERMS 


X(l) 

= 

-105.9/12.0 

X(  2  ) 

= 

-99.3/12.0 

X(  3  ) 

B 

-87.3/12.0 

X{  4  ) 

= 

-66.3/12.0 

X(  5) 

= 

72.7/12.0 

X(6) 

= 

83.2/12.0 

X(7) 

= 

91.2/12.0 

X(  8  ) 

= 

99.2/12.0 

X(  9  ) 

= 

103.2/12.0 
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BR(1)  =  0.00/12.0 

BR(2)  =  8.24/12.0 

BR(3)  =  19.76/12.0 

3R(4)  =  29.36/12.0 

BR(5)  =  31.85/12.0 

BR(6)  =  27.84/12.0 

BR{7)  =  21.44/12.0 

BR(8)  =  12.00/12.0 

BR(9)  =  0.00/12.0 

COMPUTE  AREA,  CENTROID,  AND  MOMENT  OF  INERTIA 

CALL  TRAP( 9 ,BR,X,AREA) 

DO  9  1=1,9 

VEC1(I)=X(I)*BR(I) 

VEC2(I)=X(I)*VEC1{I) 

9  CONTINUE 

CALL  TRAP( 9 ,VEC1 ,X,XAA) 

XA=XAA/AREA 

CALL  TRAP( 9 ,VEC2 ,X, lA) 

WRITE  (*,1001) 

READ  (*,*)  IRES 

OPEN( 31 ,NAME=' COEF.DAT' , STATUS= ' OLD  '  ) 

READ^31,*)  RATIO,  DELB,  XGB ,  ZGB,  XB ,  ZB 

BUOY=  WEIGHT  +  DELB 

XG=XB+XGB 

ZG=ZB+2GB 

MASS  MATRIX  COEFFICIENTS 

Bl(l,l)=  MASS  -  XUDOT 
Bl(l,2)=  0.0 
Bl(l,3)=  MASS*ZG 
Bl(l,4)=  0,0 

B1 ( 2 , 1 )=  0.0 
Bl(2,2)=  MASS  -  ZWDOT 
B1 ( 2 , 3 ) =- ( ZQDOT+MASS*XG ) 

Bl(2,4)=  0.0 

B1 ( 3 , 1 )=  MASS*ZG 

B1 ( 3 , 2 ) =- ( MWD0T+MASS*XG ) 

Bl(3,3)=  lY-MQDOT 
Bl{3,4 )=  0.0 

Bl ( 4 , 1 )=  0.0 
El{4,2)=  0.0 
Bl ( 4 , 3 )=  0.0 
Bl  ( 4  ,  4  )=  1.0 
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B2(l,l)=  IX-KPDOT 
B2 ( 1 , 2 )=  0.0 

B2 ( 1 , 3 ) =- ( KVDOT+MASS*2G ) 

B2 ( 1 , 4 )=-KRDOT 
C 

B2 ( 2 , 1 )=  0.0 
B2(2,2)=  1.0 
B2(2,3)=  0.0 
B2(2,4)=  0.0 
C 

B2 ( 3 , 1 ) =- ( YPDOT+MASS*ZG ) 

B2(3,2)=  0.0 
B2(3,3)-  MASS-YVDOT 
B2(3,4)=  MASS*XG-yRDOT 
C 

B2  (  4 , 1  )— NPDOT 
B2(4,2)=  0.0 
B2(4,3)=  MASS*XG-NVDOT 
B2(4,4)=  IZ-NRDOT 
C 

OPEN (41, NAME* ' DEOS . RES ' , STATUS* ' NEW ' ) 
OPEN(42,NAME='DEOSl.RES' , STATUS* ' NEK ' ) 

OPEN (43, NAME* ' DEOS2 . RES ' , STATUS* ' NEW ' ) 

C 

IF  (IRES.EQ.l)  GO  TO  1 
IF  (IRES.EQ.2)  GO  TO  2 
IF  (IRES.EQ.3)  GO  TO  3 
IF  (IRES.EQ.4)  GO  TO  4 

1  OPEN  ( 21, NAME* 'STl. RES ', STATUS* 'OLD' ) 

11  READ  (21,*,END=100)  DSD , THETO , UO , WO , WP 
GO  TO  5 

2  OPEN  ( 22, NAME*' ST2. RES' , STATUS* 'OLD' ) 

12  READ  ( 22 , * , END*100 )  DSD , THETO , UO , WO , KP 
GO  TO  5 

3  OPEN  ( 2 3, NAME*' ST3. RES' , STATUS* 'OLD' ) 

13  READ  ( 23 , * ,END*100 )  DSD , THETO , UO , WO , WP 
GO  TO  5 

4  OPEN  ( 24 , NAME*' ST4 .RES' , STATUS* ' OLD' ) 

14  READ  ( 24 , * ,END*100 )  DSD , THETO , UO , WO , WP 
GO  TO  5 

C 

5  THETA0*THET0*PI/180 . 0 
DS  *  DSD*PI/180.0 
DB  »  DS  *  RATIO 

DAMPING  MATRIX  COEFFICIENTS 

A1  (  1  ,  1  )*-2 . 0*U0*CD0+W0*  (XWDS^DS*XWDB’'-DB  ) 

&  +2 . 0*U0* (XDSDS*DS**2+XDBDB*DB**2 ) 

A1 ( 1 , 2 )*  2 . 0*XWW*W0+U0* ( XWDS*DS+XWDB*DB ) 
A1  (  1  ,  3  )*  (XWQ-MASS  }  *W0- (  XODS-D.'-^XODE*DB  )  *U0 
A1 ( 1 , 4 )=- (WEIGHT-BUOY )*DCOS(THErA0 ) 
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A1 ( 2 , 1 ) -  ZW*W0+2 . 0*U0* ( ZDS*DS+2DB*DB ) 

A1 (2 ,2 )«  ZW*U0-2.0*CDZ*AREA*DABS(W0) 

A1 ( 2 , 3 )-  ( ZQ+MASS )*U0+2. 0*CDZ*AREA*XA*DABS(WO ) 

A1 ( 2 , 4 ) — ( WEIGHT-BUOY ) *DSIN( THETAO ) 

C 

Al(3,l) 

A1 ( 3 , 2  ) 

A1 ( 3 , 3  ) 

& 

Al(3,4) 

& 

C 

Al(4,l)-  0.0 
Al(4,2)-  0.0 
Al(4,3)-  ■'  .0 
Al(4,4)-  0.0 
C 

A2{l,l)-=  KP*U0+(KWP-MASS*ZG)*W0 

A2 ( 1 , 2 ) =- ( ZG*WEIGHT-ZB*BUOY ) *DCOS ( THETAO ) 

A2(l,3)«  KV*U0+KVW*W0 
A2(l,4)=  (KR+MASS*ZG)*U0+KWR*W0 
C 

A2(2,l)=  1.0 
A2(2,2)=  0.0 
A2(2,3)-  0.0 
A2(2,4)=  DTAN( THETAO) 

C  • 

A2(3,l)*  YP*U0+(yWP+MASS)*W0 

A2  {  3 , 2  )*=  (WEIGHT-BUOY)  *DCOS(  THETAO  ) 

A2(3,3)=  YV*U0+YVW*W0-CDZ*AREA*DABS(W0)  4 

A2(  3,4  )«=  YWR*WO+{YR-MASS)*UO-CDZ*AREA*XA 
&  *DABS(W0) 

C 

A2(4,l)-  MASS*XG*W0+NP*U0+NWP*W0 
A2 ( 4 , 2 ) -  ( XG*WEIGHT-XB*BUOY ) *DCOS ( THETAO ) 

A2 ( 4 , 3 )-  NV*UO+NVW*WO-CDZ*AREA*XA*DABS(WO ) 

A2( 4 ,4 )-  (NR-MASS*XG)*UO+NWR*WO-CDZ*IA*DABS(WO ) 

RESTORE  B-MATRIX 

DO  71  1-1,4 
DO  72  J-1,4 

BB1(I,J)-B1(I,J) 

72  CONTINUE 

71  CONTINUE 

DO  81  1-1,4 
DO  82  J-1,4 

BB2( I , J)=B2( I , J) 

82  CONTINUE 

81  CONTINUE 


-  MW*W0+2 . 0*U0* {MDS*DS+MDB*DB ) 

-  MW*U0+2 . 0*CDZ*AREA*XA*DABS (WO ) 

-  (MQ-MASS*XG)*UO-MASS*ZG*WO 

-2.0*CDZ*1A*DABS(W0 ) 

-  ( XG*WEIGHT-XB*BUOY ) *DS1N ( THETAO ) - 
( ZG*WEIGHT-ZB*BUOY ) *DCOS ( THETAO ) 
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CALL  RGG( 4 , 4 ,Al ,BBl ,ALFRl ,ALFIl ,BETAl , 0,2ZZl , lER) 
CALL  DEGSTB(DE0S1 ,ALFR1 ,ALFI1 ,BETAl , FREQl ,WRl ,WIl ) 
C 

CALL  RGG( 4 , 4 , A2 , BB2 , ALFR2 , ALFI2 , BETA2 , 0 , ZZZ2 , lER ) 
CALL  DEGSTB ( DE0S2 , ALFR2 , ALFI2 , BETA2 , FREQ2 , WR2 , WI 2 ) 
C 

IF  (DE0S1.GE.DE0S2)  DEOS-DEOSl 
IF  ( DEOSl .LT.DEOS2 )  DEOS-DEOS2 
C 

WRITE  (41,2001)  DSD,THETO,UO,WO,WP,DEOS, 

&  DEOSl, DEOS2 

IF  (DEOS.LT.O.DO) 

&  WRITE  (42,2001)  DSD,THETO,UO,WO,WP,DEOS, 

&  DEOSl, DE052 

IF  (DEOSl .LT. 0 .DO ) 

&  WRITE  (43,2001)  DSD , THETO , UO , WO , WP , DEOS , 

&  DEOSl, DEOS2 

C 

IF  ( IRES.EQ.l )  GO  TO  11 

IF  (IRES.EQ.2)  GO  TO  12 

IF  (IRES.EQ.3)  GO  TO  13 

IF  (IRES.EQ.4)  GO  TO  14 

C 

100  STOP 
1001  FORMAT 
& 

2001  FORMAT 

2002  FORMAT 
END 

C 

SUBROUTINE  DEGSTB ( DEOS ,ALFR , ALFI , BETA , OMEGA , WR ,WI ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  ALFR( 4 ) ,ALFI ( 4 ) ,BETA( 4 ) ,WR( 4 ) ,WI ( 4 ) 

DO  1  1=1,4 

WR( I )=ALFR( I )/BETA( I ) 

WI ( I )=ALFI ( I )/BETA( I ) 

1  CONTINUE 
DEOS=-l . OE+10 
DO  2  1=1,4 

IF  (WR( I ) .LT.DEOS )  GO  TO  2 
DEOS=WR( I ) 

IJ-I 

2  CONTINUE 
OMEGA=WI ( I J ) 

OMEGA=DABS ( OMEGA ) 

RETURN 
END 

C 

SUBROUTINE  TRAP ( N , A , B , OUT ) 


(  '  ENTER  THE  RESPONSE  DATA  FILE  DESIRED 
(1,2,3,  OR  4 )  ' ) 

(8E15.5) 

(F10.3) 
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NUMERICAL  INTEGRATION  ROUTINE  USING 
THE  TRAPEZOIDAL  RULE 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  A( 1 ) ,B( 1 ) 

N1«N-1 
OUT=0 . 0 
DO  1  1=1, N1 

OUT1=0.5* (A( I )+A( I+l ) )* (B{ I+l )-B(I ) ) 
OUT  =OUT+OUTl 
1  CONTINUE 
RETURN 
END 


a 


a 
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